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Preface 



Students of science and engineering often encounter a gap between their 
calculus text and the applied mathematical books they need to read soon 
after. What is that gap? A calculus student might have learned some 
calculus without learning how to read mathematics. A science major might 
be faced with an applied text which assumes the ability to read a somewhat 
encyclopedic mathematics and its application. I often find myself teaching 
students who are somewhere between these two states. 

These notes are mostly from lectures I have given to bridge that gap, while 
teaching the Engineering Mathematics courses at Cornell University. The 
notes could be used for an introductory unified course on ordinary and par- 
tial differential equations. However, this is not a textbook in the usual sense, 
and certainly not a reference. The stress is on clarity, not completeness. It 
is an introduction and an invitation. Most ideas are taught by example. 
The notes have been available for many years from my web page. 

For fuller coverage, see any of the excellent books: 

• Agnew, Ralph Palmer, Differential Equations, McGraw-Hill, 1960 

• Churchill, Ruel V., Fourier Series and Boundary Value Problems, Mc- 
Graw Hill, 1941 

• Hubbard, John H., and West, Beverly H., Differential Equations, a 
Dynamical Systems Approach, Parts 1 and 2, Springer, 1995 and 1996 

• Lax, P., Burstein, S., and Lax, A., Calculus with Applications and 
Computing: Vol 1, Springer, 1983 

• Seeley, Robert T., An Introduction to Fourier Series and Integrals, 
W.A.Benjamin, 1966 

Some of the exercises have the format "What's rong with this?". These 
are either questions asked by students or errors taken from test papers of 
students in this class, so it could be quite beneficial to study them. 

Robert E. Terrell 

January 1997 September 2005 April 2009 
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1 Introduction to Ordinary Differential Equations 



Today: An example involving your bank account, and nice pictures called 
slope fields. How to read a differential equation. 

Welcome to the world of differential equations! They describe many pro- 
cesses in the world around you, but of course we'll have to convince you of 
that. Today we are going to give an example, and find out what it means 
to read a differential equation. 



1.1 The Banker's Equation 

A differential equation is an equation which contains an unknown function 
and one of its unknown derivatives. Here is a differential equation. 

d J = -02% 
dt y 

It doesn't look too exciting does it? Really it is, though. It might for 
example represent your bank account, where the balance is y at a time t 
years after you open the account, and the account is earning 2.8% interest. 
Regardless of the specific interpretation, let's see what the equation says. 
Since we see the term dy/dt we can tell that y is a function of t, and that 
the rate of change is a multiple, namely .028, of the value of y itself. We 
definitely should always write y(t) instead of just y, and we will sometimes, 
but it is traditional to be sloppy. 

For example, if y happens to be 2000 at a particular time t, the rate of 
change of y is then .028(2000) = 56, and the units of this rate in the bank 
account case are dollars/year. Thus y is increasing, whenever y is positive. 

Practice: What do you estimate the balance will be, roughly, a year from 
now, if it is 2000 and is growing at 56 dollars/year? 

This is not supposed to be a hard question. By the way, when I ask a 
question, don't cheat yourself by ignoring it. Think about it, and future 
things will be easier. I promise. 

Later when y is, say, 5875.33, its rate of change will be .028(5875.33) = 
164.509 which is much faster. We'll sometimes refer to y' = .028y as the 
banker's equation. 
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Do you begin to see how you can get useful information from a differential 
equation fairly easily, by just reading it carefully? One of the most important 
skills to learn about differential equations is how to read them. For example 
in the equation 

y' = .02% - 10 

there is a new negative influence on the rate of change, due to the —10. This 
— 10 could represent withdrawals from the account. 

Practice: What must be the units of the -10? 

Whether the resulting value of y' is actually negative depends on the current 
value of y. That is an example of "reading" a differential equation. As a 
result of this reading skill, you can perhaps recognize that the banker's 
equation is very idealized: It does not account for deposits or changes in 
interest rate. It didn't account for withdrawals until we appended the —10, 
but even that is an unrealistic continuous rate of withdrawal. You may wish 
to think about how the equation might be modified to include those things 
more realistically. 

1.2 Slope Fields 

It is significant that you can make graphs of the solutions sometimes. In 
the bank account problem we have already noticed that the larger y is, the 
greater the rate of increase. This can be displayed by sketching a "slope 
field" as in Figure 1.1. Slope fields are done as follows. First, the general 
form of a differential equation is 

y' = f(y,t) 

where y(t) is the unknown function, and / is given. To make a slope field 
for this equation, choose some points (y, t) and evaluate / there. According 
to the differential equation, these numbers must be equal to the derivative 
y', which is the slope for the graph of the solution. These resulting values 
of y' are then plotted using small line segments to indicate the slopes. For 
example, at the point t = 6, y = 20, the equation y' = .02% says that the 
slope must be .028(20) = .56. So we go to this point on the graph and 
place a mark having this slope. Solution curves then must be tangent to the 
slope marks. This can be done by hand or computer, without solving the 
differential equation. 
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Figure 1.1 The slope field for y = .028y as made by df ield. In Lecture 7 
we will discuss where to get df ield. Some solutions are also shown. 

Note that we have included the cases y = and y < in the slope field even 
though they might not apply to your bank account. (Let us hope.) 

Practice: Try making a slope field for y' — y + t. 

There is also a way to explicitly solve the banker's equation. Assume we are 
looking for a positive solution. Then it is alright to divide the equation by 
y, getting 




= .028 



Then integrate, using the chain rule: 

ln(y) = .028t + c 

where c is some constant. Then 

y (t) = e- 028t+c = Cl e- 02St 

Here we have used a property of the exponential function, that e a+b = e a e b , 
and set c\ = e c . The potential answer which we have found must now be 
checked by substituting it into the differential equation to make sure it really 
works. You ought to do this. Now. You will notice that the constant c\ can 
in fact be any constant, in spite of the fact that our derivation of it seemed 
to suggest that it be positive. 

This is how it works with differential equations: It doesn't matter how you 
solve them, what is important is that you check to see whether you are right. 



Some kinds of answers in math and science can't be checked very well, but 
these can, so do it. Even guessing answers is a highly respected method! if 
you check them. 

Your healthy scepticism may be complaining to you at this point about any 
bank account which grows exponentially. If so, good. It is clearly impossible 
for anything to grow exponentially forever. Perhaps it is reasonable for a 
limited time. The hope of applied mathematics is that our models will be 
idealistic enough to solve while being realistic enough to be worthwhile. 

The last point we want to make about this example concerns the constant 
ci further. What is the amount of money you originally deposited, y(0)? Do 
you see that it is the same as ci? That is because y(0) = ere and e° = 1. If 
your original deposit was 300 dollars then c\ = 300. This value y(0) is called 
an "initial condition", and serves to pick the solution we are interested in 
out from among all those which might be drawn in Figure 1.1. 

Example: Be sure you can do the following kind 

x = —3a; 
x(0) = 5 

Like before, we get a solution x(t) = ce~ 3 *. Then a;(0) = 5 = ce° = c so 
c = 5 and the answer is 

x(t) = 5e~ 3t 

Check it to be absolutely sure. 
Problems 

1. Make slope fields for x' — x, x' — —x, x' = x + t, x' = —x + cost, x' = t. 

2. Sketch some solution curves onto the given slope field in Figure 1.2 below. 

3. What general fact do you know from calculus about the graph of a function y if y' > 0? 
Apply this fact to any solution of y' = y — y 3 : consider cases where the values of y lie in 
each of the intervals (-co, — 1),(— 1, 0), (0, 1), and (1, oo). For each interval, state whether 
y is increasing or decreasing. 

4. (continuing 3.) If y' = y — y 3 and y (0) = | , what do you think will be the lim t ^oo y(t)7 
Make a slope field if you're not sure. 

5. Reconsider the banker's equation y' = .028y. If the interest rate is 3% at the beginning 
of the year and expected to rise linearly to 4% over the next two years, what would you 
replace .028 with in the equation? You are not asked to solve the equation. 

6. In y' = .028y — 10, suppose the withdrawals are changed from $10 per year continuously, 
to $200 every other week. What could you replace the —10 by? This is a difficult rhetorical 
question only- something to think about. Do you think it would be alright to use a smooth 
function of the form a cos bt to approximate the withdrawals? 



4 



x' = x A 2-sin(t) 
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Figure 1.2 The slope field for Problem 2, as drawn by df ield. 



2 A Gallery of Differential Equations 

Today: A gallery is a place to look and to get ideas, and to find out how 
other people view things. 

Here is a list of differential equations, as a preview of things to come. Unlike 
the banker's equation y' = 0.028y, not all differential equations are about 
money. 

1. y' = .028(600 — y) This equation is a model for the heating of a pizza 
in a 600 degree oven. Of course the 0.028 is there just for comparison with 
the banker's interest. The physical law involved is called, believe it or not, 
"Newton's Law of Cooling". We'll do it in Lecture 8. 

2. Newton had other laws as well, one of them being the ll F = ma" law 
of inertia. You might have seen this in a physics class, but not realised 
that it is a differential equation. That is because it concerns the unknown 
position of a mass, and the second derivative a, of that position. In fact, 
the original differential equation, the very first one over 300 years ago, was 
made by Newton for the case in which F was the gravity force between 
the earth and moon, and m was the moon's mass. Before that, nobody 
knew what a differential equation was, and nobody knew that gravity had 
anything to do with the motion of the moon. They thought gravity was 
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what made their physics books heavy. You have probably heard a garbled 
version of the story in which Newton saw or was hit by a falling apple while 
thinking about gravity What is frequently omitted is that he looked up to 
see whence it had fallen, and saw the moon up there, behind the branches. 
Thus it occured to him that there might be a connection. Maybe his thought 
process was something like, "gravity, gravity, gravity, apple, gravity, gravity, 
apple, moon, gravity, apple, moon, gravity, moon, gravity, moon, gravity. 
Oh! Moon Gravity!" 

There is a joke by the famous physicist Richard Feynmann, who said that 
once people thought there were angels up there pushing against the side of 
the moon to keep it moving around the earth. This is of course ridiculous, 
he said. After Newton everyone knows that the angels are in fact on the far 
side of the moon pushing it toward us! 

Anyway Newton's law often leads to second order differential equations. 

,2 

Practice: The simplest second order equations look like = —x. You 
ought to be able to guess some solutions to this equation. What are they? 

3. There are also Partial Differential Equations. This does not mean that 
somebody forgot to write out the whole thing! It means that the unknown 
functions depend on more than one variable, so that partial derivatives show 
up in the equation. One example is ut = u xx where the subscripts denote 
partial derivatives. This is called the heat equation, and u(x, t) is the tem- 
perature at position x and time t, when heat is allowed to conduct only 
along the x axis, as through a wall or along a metal bar. It concerns a dif- 
ferent aspect of heat than does Newton's law of cooling, and we will discuss 
it more in Lectures 11 and 27. 

4. uu = u xx looks sort of like the heat equation, but is very different be- 
cause of the second time derivative. This is the wave equation, which is 
about electromagnetic waves (wireless), music, and water waves, in decreas- 
ing order of accuracy. The equation for vibrations of a drum head is the two 
dimensional wave equation, uu = u xx + u yy . In Lecture 34 we will derive a 
special polar coordinate version of that, uu = u rr + \u r . We'll use that to 
describe some of the sounds of a drum. 

5. There are others which we won't study but some of the ideas we use can 
be applied to them. The reason for mentioning these here is to convince 
you that the earlier statement about "many processes" in the world around 
you is correct. If you bang your fist on the table top, and the table top is 
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somewhat rigid, not like a drum head, then the sound which comes out is 
caused by vibrations of the wood, and these are described by solutions of 
u u = -{u 

xxxx H~ ^Uxxyy ~\~ Hyyyy) * H6F6 u(x, y, t) is the vertical deflection of 
the table top while it bends and vibrates. Can you imagine that happening 
at a small scale? We will in fact do a related equation u tt = — u xxxx for 
vibrations of a beam in Lecture 31. 

At a much smaller scale, the behavior of electrons in an atom is described 
by Schrodinger's equation iiit = —{u xx + u yy + u zz ) +V(x,y, z)u. In this 
case, u(x, y, z, t) is related to the probability that the electron is at (x, y, z) 
at time t. At the other end of the scale there is an equation we won't 
write down, but which was worked out by Einstein. Not E = mc 2 , but a 
differential equation. That must be why Einstein is so famous. He wrote 
down a differential equation for the whole universe. 

Problems 

1. Newton's gravity law says that the force between a big mass at the origin of the x axis 
and a small mass at point x(t) is proportional to x~ 2 . How would you write the F — ma 
law for that as a differential equation? 

2. We don't have much experience with the fourth derivatives mentioned above. Let 
u(x,t) — t 2 x 2 — x 3 + x 4 . Is u XX xx positive or negative? Does this depend on the value of 

x? 

3. In case you didn't do the Practice item, what functions do you know about from 
calculus, that are equal to their second derivative? the negative of their second derivative? 



3 Introduction to Partial Differential Equations 



Today: A first order partial differential equation. 



Here is a partial differential equation, sometimes called a transport equation, 
and sometimes called a wave equation. 

dw g dw g 
dt dx 



Practice: We remind you that partial derivatives are the rates of change 
holding all but one variable fixed. For example 

(x - y 2 x + 2yt) = 2y, ^- {x - y 2 x + 2yt) = 1 - y 2 
What is the y partial? 
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Our PDE is abbreviated 



w t + 3w x = 







You can tell by the notation that w is to be interpreted as a function of 
both t and x. You can't tell what the equation is about. We will see that 
it can describe certain types of waves. There are water waves, electromag- 
netic waves, the wavelike motion of musical instrument strings, the invisible 
pressure waves of sound, the waveforms of alternating electric current, and 
others. This equation is a simple model. 

Practice: You know from calculus that increasing functions have positive 
derivatives. In Figure 3.1 a wave shape is indicated as a function of x at 
one particular time t. Focus on the steepest part of the wave. Is w x positive 
there, or negative? Next, look at the transport equation. Is Wt positive 
there, or negative? Which way will the steep profile move next? 
Remember how important it is to read a differential equation. 

3.1 A Conservation Law 

We'll derive the equation as one model for conservation of mass. You might 
feel that the derivation of the equation is harder than the solving of the 
equation. 

We imagine that w represents the height of a sand dune which moves by the 
wind, along the x direction. The assumption is that the sand blows along 
the surface, crossing position (x,w(x, t)) at a rate proportional to w. The 
proportionality factor is taken to be 3, which has dimensions of velocity, like 
the wind. 



Figure 3.1 The wind blows sand along the surface. Some enters the segment 
(a, a+h) from the left, and some leaves at the right. The net difference causes 
changes in the height of the dune there. 

The law of conservation of sand says that over each segment (a, a + h) you 





a 



a+h 



have 




a+h 
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That is the time rate of the total sand on the left side, and the sand flux on 
the right side. Divide by h and take the limit. 

Practice: 1. Why is there a minus sign on the right hand term? 

2. What do you know from the Fundamental Theorem of Calculus about 



The limit we need is the case in which / is w t (x,t). 

We find that wt(a, t) = —3w x (a, t). Of course a is arbitrary. That concludes 
the derivation. 

3.2 Traveling Waves 

When you first encounter PDE, it can appear, because of having more than 
one independent variable, that there is no reasonable place to start working. 
Do I try t first, x, or what? In this section we'll just explore a little. If we 
try something that doesn't help, then we try something else. 

Practice: Find all solutions to our transport equation of the form 

w(x, t) = ax + bt 

In case that is not clear, it does not mean 'derive ax + bt somehow'. It means 
substitute the hypothetical w(x,t) = ax + bt into the PDE and see whether 
there are any such solutions. What is required of a and b? 

Those practice solutions don't look much like waves. 



Figure 3.1 Graphs of cos(a;), cos(x — 1.5), and cos(a; — 3). If you think of 
these as photographs superimposed at three different times, then you can 
see that it moves. What values of t do you associate with these graphs, and 
how fast does the wave move? 

Lets try something more wavey. 
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Practice: Find all solutions to our wave equation of the form w(x, t) = 
ccos(aa: + bt) 

So far, we have seen a lot of solutions to our transport equation. Here are a 
few of them: 

w(x, t) = x — 3t 

w(x,t) = -2.1x + 6M 

w(x,i) = 40cos(5x - 15i) 
2 

w(x,t) = —- cos(8x - 24i) 

For comparison, that is a lot more variety than we found for the banker's 
equation. Remember that the only solutions to the ODE y' = .028y are 
constant multiples of e - 028 *. Now lets go out on a limb. Since our wave 
equation allows straight lines of all different slopes and cosines of all different 
frequencies and amplitudes, maybe it also allows other things too. 

Try 

w(x, t) = f(x — 3t) 

where we won't specify the function / yet. Without specifying / any further, 
we can't find the derivatives we need in any literal sense, but can apply the 
chain rule anyway. The intention here is that / ought to be a function of 
one variable, say s, and that the number x — 3t is being inserted for that 
variable, s = x — 3t. The partial derivatives are computed using the chain 
rule, because we are composing / with the function x — 3t of two variables. 
The chain rule here looks like this: 

dw dfds , 

Wt= m = dSdi = - 3f{x - 3t) 

Practice: Figure out why w x = f'(x - St). 

Setting those into the transport equation we get 

w t + 3w x = -3f + 3f = 

That is interesting. It means that any differentiable function / gives us a 
solution. Any dune shape is allowed. You see, it doesn't matter at all what 
/ is, as long as it is some differentiable function. 
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Don't forget: differential equations are a model of the world. They are not 
the world itself. Real dunes cannot have just any shape / whatsoever. They 
are more specialized than our model. 

Practice: Check the case f(s) = 22sin(s) — 10sin(3s). That is, verify 
that 

w(x, t) = 22 sin(a; - 3t) - 10 sin(3:r - 9t) 
is a solution to our wave equation. 

Problems 

1. Work all the Practice items in this lecture if you have not done so yet. 

2. Find a lot of solutions to the wave equation 

yt -5y x = 

and tell which direction the waves move, and how fast. 

3. Check that w(x, t) = 1/(1 + (x — 3i) 2 ) is one solution to the equation w t + 3w x = 0. 

4. What does the initial value w(x, 0) look like in problem 3, if you graph it as a function 
of x? 

5. Sketch the profile of the dune shapes w(x, 1) and w(x, 2) in problem 3. What is 
happening? Which way is the wind blowing? What is the velocity of the dune? Can you 
tell the velocity of the wind? 

6. Solve u t + u x — if we also want to have the initial condition u(x, 0) = | cos(2a;). 

7. As in problem 6, but with initial condition u(x, 0) = | cos(2:e) + | sin(4a;). Sketch the 
wave shape for several times. 
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4 The Logistic Population Model 



Today: The logistic equation is an improved model for population growth. 

We have seen that the banker's equation y' = .028y has exponentially grow- 
ing solutions. It also has a completely different interpretation from the bank 
account idea. Suppose that you have a population containing about y(t) in- 
dividuals. The word "about" is used because if y = 32.51 then we will have 
to interpret how many individuals that is. Also the units could be, say, 
thousands of individuals, rather than just plain individuals. The population 
could be anything from people on earth, to deer in a certain forest, to bacte- 
ria in a certain Petrie dish. We can read this differential equation to say that 
the rate of change of the population is proportional to the number present. 
That perhaps captures some element of truth, yet we see right away that no 
population can grow exponentially forever, since sooner or later there will 
be a limit imposed by space, or food, or energy, or something. 

The Logistic Equation 

Here is a modification to the banker's equation that overcomes the previous 
objection. 

| = . 028,(1 -y) 

In order to understand why this avoids the exponential growth problem we 
must read the differential equation carefully. Remember that I said this is 
an important skill. 

Here we go. You may rewrite the right-hand side as .028(y — y 2 )- You know 
that when y is small, y 2 is very small. Consequently the rate of change 
is still about .028y when y is small, and you will get exponential growth, 
approximately. After this goes on for a while, it is plausible that the y 2 
term will become important. In fact as y increases toward 1 (one thousand 
or whatever), the rate of change approaches 0. 

For simplicity we now dispense with the .028, and for flexibility introduce a 
parameter a, and consider the logistic equation y' = y(a — y). Therefore if 
we make a slope field for this equation we see something like Figure 4.1. 
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Figure 4.1 A slope field for the Logistic Equation. Note that solutions 
starting near have about the same shape as exponentials until they get 
near a. (figure made using df ield) 

The solutions which begin with initial conditions between and a evidently 
grow toward a as a limit. This in fact can be verified by finding an explicit 
formula for the solution. Proceeding much as we did for the bank account 
problem, first write 

k dy = dt 

ay - y z 

To make this easier to integrate, we'll use a trick which was discovered by a 
student in this class, and multiply first by y~ 2 jy~ 2 . Then integrate 

\n(ay~ l - l)=t + c 

a 

The integral can be done without the trick, using partial fractions, but that 
is longer. Now solve for y(t) 

ay- 1 - 1 = e - a (' +c ) = cie~ at 

i \ a 

1 + c\e at 

These manipulations would be somewhat scary if we had not specified that 
we are interested in y values between and a. For example, the In of a 
negative number is not defined. However, we emphasize that the main point 
is to check any formulas found by such manipulations. So let's check it: 

a 2 ae- at 

y 



[1 + cie~ at ) 2 
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We must compare this expression to 



a 



2 



a 



2 



2 ( 1+Cie -at)_l 

(l + cie-«*) 2 



ay -y 2 = 



i + ae~ at (l + ae- at ) 2 



= a' 



You can see that this matches y'. Note that the value of c\ is not restricted 
to be positive, even though the derivation above may have required it. We 
have seen this kind of thing before, so checking is very important. The only 
restriction here occurs when the denominator of y is 0, which can occur if 
ci is negative. If you stare long enough at y you will see that this does not 
happen if the initial condition is between and a, and that it restricts the 
domain of definition of y if the initial conditions are outside of this interval. 
All this fits very well with the slope field above. In fact, there is only one 
solution to the equation which is not contained in our formula. 

Practice: Can you see what it is? 



1. Suppose that we have a solution y(t) for the logistic equation y' — y(a — y). Choose 
some time delay, say 3 time units to be specific, and set z(t) = y(t — 3). Is z(t) also a 
solution to the logistic equation? 

2. The three 'S'-shaped solution curves in Figure 4.1 all appear to be exactly the same 
shape. In view of Problem 1, are they? 

3. Prof. Verhulst made the logistic model in the mid-1800s. The US census data from the 
years 1800, 1820, and 1840, show populations of about 5.3, 9.6, and 17 million. We'll need 
to choose some time scale ti in our solution y(t) = a(l + cie _ai ) _1 so that t = means 
1800, t — t\ means 1820, and t = 2ti means 1840. Figure out ci, ti, and a to match the 
historical data, warning: The arithmetic is very long. It helps if you use the fact that 
e -o-2t! = ( e -oti)2_ ANSWER: ci = 36.2, ti = .0031, and a = 197. 

4. Using the result of Problem 3, what population do you predict for the year 1920? The 
actual population in 1920 was 106 million. The Professor was pretty close wasn't he? He 
was probably surprised to predict a whole century. 

Example: Census data for 1810, 1820, and 1830 show populations of 7.2, 
9.6, and 12.8 million. Trying those, it turned out that I couldn't fit the 
numbers due to the numerical coincidence that (7.2)(12.8) = (9.6) 2 . That is 
where I switched to the years in Problem 3. This shows that the fitting of 
real data to a model is nontrivial. Mathematicians like the word 'nontrivial'. 



Problems 



5 Separable Equations 



Today: Various useful examples. 
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Equations of the form 

include the banker's and logistic equations and some other useful equations. 
You may attack these by "separating variables": Write -J^j = g(t)dt, then 
try to integrate and solve for x{t). Remember that this is what we did with 
the banker's and logistic equations. 

Example: 

dx 

■dt =tX 

= tdt, In \x\ = y + c, x(t) = cie5* 2 , check : x' = yCieJ* 2 = tx 

The skeptical reader will wonder why I said "try" to integrate and solve for 
x. The next two examples illustrate that there is no guarantee that either 
of these steps may be completed. 

Example: 

dx _ 1 
dt ~ l + 3x 2 

(1 + 3a; 2 ) dx — dt,x + x 3 — t + c. Here you can do the integrals, but you 
can't very easily solve for x. 



Example: 

— - e-* 2 
dt 

_ ,2 

dx = e dt. Here you can't do the integral. 
The Leaking Bucket 

(This example comes from the Hubbard and West book. See preface.) 

We're about to write out a differential equation to model the depth of water 
in a leaking bucket, based on plausible assumptions. We'll then try to solve 
for the water depth at various times. In particular we will try to reconstruct 
the history of an empty bucket. This problem is solvable by separation of 
variables, but it has far greater significance that the problem has more than 
one solution. Physically if you see an empty bucket with a hole at the 
bottom, you just can't tell when it last held water. 

Here is a crude derivation of our equation which is not intended to be physi- 
cally rigorous. We assume that the water has depth y{t) and that the speed 
of a molecule pouring out at the bottom is determined by conservation of 
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energy. Namely, as it fell from the top surface it gained potential energy 
proportional to y and this became kinetic energy at the bottom, propor- 
tional to (y 1 ) 2 . Actually y' is the speed of the top surface of the water, not 
the bottom, but these are proportional in a cylindrical tank. Also we could 
assume a fraction of the energy is lost in friction. In any event we get (y') 2 
proportional to y, or 

y' = -ay/y 

The constant is negative because the water depth is decreasing. We'll set 
a = 1 for convenience. 

The slope field looks like the figure. 
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Figure 5.1 Depth of water in the leaking bucket. 



Separation of variables gives, if y > 0, 

-y- 1 ' 2 dy = dt 
2y 1 ' 2 = -t + c 

so here we must have t < c, otherwise the square root would be negative. 



y 



Note that if t = c then y = so the derivation is no good there. Comparing 
our formula to the slope field we see what must happen: 



y(t) 



(^) 2 ift<c 
otherwise 
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The final point here is that if our initial condition is, say, y(0) = 0, then any 
choice of c < gives an acceptable solution. This means that if the bucket 
is empty now, it could have become empty at any time in the past. Again, 
we consider this kind of behavior to be unusual, in the sense that we don't 
usually notice or want examples of nonunique behavior. Curiously in this 
problem, nonuniqueness fits reality perfectly. 

Example: For 

^ = ^x 3 with x(0) = -3 

you calculate x~ 3 dx = \ dt, — \x~ 2 = \t + c, x~' 2 — — t — 2c, (so t < —2c). 
Then at t = we get (-3)~ 2 = -0 - 2c, 2c = -1/9, x~ 2 = -t + 1/9, 
x = ±(— t + 1/9) -1 / 2 . Now recheck the initial conditions to determine the 
sign: —3 = ±(0 + 1/9)~ 1//2 , therefore use the minus sign: 

x = -(-t + l/9r 1/2 for (it < 1/9). 

Check it: -3 = -(1/9)" 1/2 is ok. x' = -(-§)(-*+ f) - * = \x 3 is also ok. 

Problems 

1. Repeat the previous example, x' — \x 3 , for the case x(0) = 3. 

2. Solve x' — — \x with initial condition a;(0) = —3. 

3. Solve x' — — |x 3 with initial condition x(0) — 3. 

4. Sketch slope fields for problems 1 and 3. 

5. Solve y' = y — y . 

6. Solve x + t 2 x = 0. 

7. Solve x + t 3 x = 0. 

8. Solve x = t 4 x. 

9. Solve x' = .028(1 + cost)x. 

10. Solve x' = a(t)x. 

11. What's rong with this? x' = x 2 , j§ = dt, ln(a; 2 ) = t, x 2 = e* + c. There are at least 
two errors. 

6 Separable Partial Differential Equations 

Today: Another use of the term "separable" . This time: PDE. You could 
wait to read this until after you have read Lecture 28, if you plan to focus 
on ODE. 
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Traditionally a distinction in terminology is made between ordinary DE with 
the variables separable, and the PDE solvable by separation of variables. 
That sounds confusing, doesn't it? It is not just the sound. The method 
looks different too. Here is an example using a transport equation 

ut + cu x = 

The method of separation of variables goes like this. Suppose we look for 
special solutions of the form 

u(x,t) = X(x)T(t) 

This does not mean that we think all solutions would be like this. We hope 
to find that at least some of them are like this. Substitute this into our PDE 
and you get 

X(x)T'(t) + cX'(x)T(t) = 

Now what? The problem is to untangle the x and t somehow. Notice that if 
you divide this equation by X(x) then the first term will at least no longer 
contain x. Lets do that, assuming temporarily that we are not dividing by 
0. 

X'(x)T(t) 

Now what? If you divide by T(t) you can clean up the second term without 
hurting the first one. Lets do that too. 

T(t) +C X{x) 

The first term only depends on t, but the second term does not depend on 
t at all. The more you think about that, the stranger it becomes. We want 
to argue that this forces the two terms to be constant. Write 

T'{t) _ _ X'(x) 

m~~ c xjx) 

For example, when x = the right hand side has a certain value, perhaps 
x' (o) 

it is 5.37. Lets write xTo/ = r wnere —cr is 5.37 or whatever. So the left 
side has to be 5.37 too. Next change x to something else, say x = 1.5. This 
doesn't affect the left side, 5.37, because there is no x in the left side. But 
then the right side has to still be 5.37. You see, the equation means that 
each side is the same constant. Calling that constant — cr, we have that 

T'(t) = -crT(t) and X'(x) = rX{x) 
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where r is some number. Thus we have split our PDE into two ODEs, 
at least for the purpose of finding a few solutions, maybe not all of them. 
Fortunately after all these new ideas, these ODEs are equations that we 
already know how to solve. We find T(t) = e~ crt and X(x) = e rx or any 
constant multiples of these. 

We have tentatively found a lot of solutions of the form 

u(x, t) = e r{x ~ ct) 

and constant multiples of these for our PDE. Among these are e x ~ ct , e ct ~ x , 
16e 25ct " 25:r , etc. 

It is always important to check any proposed solutions, particularly since 
our derivation of these was so new a method. In fact, going a little further, 
lets try a linear combination of two of them 

u(x, t) = e x ~ ct + lQe 25ct ~ 25x 

We get u t = -ce x ~ ct + 25c • l6e 25c *- 25:r and u x = e x ~ ct - 25 • I6e 25c *- 25:E . So 
Ut + cu x = right enough. Apparently you can add any number of these 
things and still have a solution. 

That worked really well in a certain sense. But on the other hand, those are 
awfully big dunes! because of the exponentials. 

If you have read about traveling waves on page 8 you know a lot of other 
solutions to the transport equation. So we won't worry about the fact that 
the separation argument only gave a few special solutions. Our purpose was 
to introduce the method. 

Two Kinds of Separating 

Finally lets compare the 'separable variables' in ODE and 'separation of 
variables' in PDE. This is rather strange, but worth thinking about. When 
we separate the ODE y' + cy = we write 

d y a + 

— = —cat 

y 

and we still consider the two sides to be dependent on each other, because 
we are looking for j/asa function of t. This has to be contrasted with our 

T X' 
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above, for the conservation law case. Here we view x and t as independent 
variables, and T and X as unrelated functions, so the approach is different. 

But the strangeness of this comparison emphasizes the fact that our deriva- 
tions are a kind of exploration, not intended as proof of anything. The plan 
is: we explore to find the form of solutions, but then we check them to be 
sure. That way, anything bogus won't botherus. 

Problems 

1. Above when we found the ODEs T'(t) = —crT(t) and X'(x) = rX(x) we wrote down 
some answers T(t) = e~ crt and X(x) = e rx pretty fast. If you aren't sure about those, 
solve these for practice. 

• t=y 

• £ = -32/ 

• T'(t) + 5T(t) = 

• /'( s ) = c /( s ) where c is a constant. 

2. Try to find some solutions to the equation u t + u x + u — by the method of separation 
of variables. Do any of your solutions have the form of traveling waves? Do they all travel 
at the same speed? 

3. The equation u t + cu x + u 2 = can't be solved by separation of variables. Try it, and 
explain why that doesn't work. 

4. In problem 3 however, there are some traveling wave solutions. Find them, if you want 
a challenge. But before you start, observe carefully that there is one special wave speed 
which is not possible. Reading the differential equation carefully, can you tell what it is? 
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7 Existence and Uniqueness and Software 



Today: We learn that some equations have unique solutions, some have 
too many, and some have none. Also an introduction to some of the available 
software. 

We have seen an initial value problem, the leakey bucket problem, for which 
several solutions existed. This phenomenon is not generally desirable, be- 
cause if you are running an experiment you would like to think that the same 
results will follow from the same initial conditions each time you repeat the 
experiment. 

On the other hand, it can also happen that an innocent-looking differential 
equation has no solution at all, or a solution which is not defined for all 
time. As one example, consider 

{x') 2 + x 2 + 1 = 

This equation has no real solution. See why? For another example consider 

x' = x 2 

By separation of variables we can find some solutions, such as x(t) = 
but this function is not defined when t = 3. In fact we should point out 
that this formula 3^ has to be considered to define two functions, not 
one, the domains being (—00, 3) and (3, 00) respectively. The reason for 
this distinction is simply that the solution of a differential equation has by 
definition to be differentiable, and therefore continuous. 

We therefore are interested in the following general statement about what 
sorts of equations have solutions, and when they are unique, and how long, 
in time, these solutions last. It is called the "Fundamental Existence and 
Uniqueness Theorem for Ordinary Differential Equations", which big title 
means, need I say, that it is considered to be important. 

7.1 The Existence Theorem 

The statement of this theorem includes the notion of partial derivative, 
which some may need practice with. So we remind you that the partial 
derivative of a function of several variables is defined to mean that the 
derivative is constructed by holding all other variables constant. For exam- 
ple, if f(t, x) = x 2 t — cos(i) then || = 2xt and ^ = x 2 + sin(t). 
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Existence and Uniqueness Theorem Consider an initial 
value problem of the form 

x'(t) = f(x,t) 

x(to) = x 

where /, to, and xo are given. Suppose it is true that / and S| are 
continuous functions of t and x in at least some small rectangle 
containing the initial condition (xo,io)- Then the conclusion is 
that there is a solution to the problem, it is defined at least for a 
small amount of time both before and after to, and there is only 
one such solution. There is no guarantee about how long a time 
interval the solution is defined for. 

Now let's see what went wrong with some of our previous examples. For the 
leakey bucket equation, we have f(t,x) = —\fx and ^ = — These are 
just fine and dandy as long as xo > 0, but when xo = you can't draw a 
rectangle around the initial conditions in the domain of /, and even worse, 
the partial derivative is not defined because you would have to divide by 0. 
The conclusion is that the Theorem does not apply. 

Note carefully what that conclusion was! It was not that existence or unique- 
ness are impossible, but only that the theorem doesn't have anything to say 
about it. When this happens you have to make a more detailed analysis. A 
theorem is really like the guarantee on your new car. It is guaranteed for 
10,000 miles, but hopefully it does not just collapse and fall to pieces all 
over the road as soon as you reach the magic number. Similarly some of the 
conclusions of a theorem can still be true in a particular example even if the 
hypotheses do not hold. To see it you can't apply the theorem but you can 
make a detailed analysis of your example. So here in the case of the leakey 
bucket we do have existence anyway. 

For the first example above, (x') 2 + x 2 + 1 = 0, the equation is not at all 
in the form required for the theorem. You could perhaps rearrange it as 
x' = \J — 1 — x 2 but this hurts my eyes-the right side is not even defined 
within the realm of real numbers in which we are working. So the theorem 
again says nothing about it. 

For the other example above, x' = x 2 , the form is alright and the right side is 
fit, x) = x 2 . This is continuous, and ^| = 2x, which is also continuous. The 
theorem applies here, but look at the conclusion. There is a solution defined 
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for some interval of t around the initial time, but there is no statement about 
how long that interval is. So our original puzzlement over the short duration 
of the solution remains. If you were a scientist working on something which 
might blow up, you would be glad to be able to predict when or if the 
explosion will occur. But this requires a more detailed analysis in each 
case-there is no general theorem about it. 

A very important point about uniqueness is the implication it has for our 
pictures of slope fields and solution curves. Two solution curves cannot cross 
or touch, when uniqueness applies. 

Practice: Look again at Figure 5.1. Can you see a point at which two 
solution curves run into each other? Take that point to be {xo,to). Read 
the Fundamental Theorem again. What goes wrong? 

7.2 Software 

In spite of examples we have seen so far, it turns out that it is not possible 
to write down solution formulas for most differential equations. This means 
that we have to draw slope fields or go to the computer for approximate 
solutions. We soon will study how approximate solutions can be computed. 
Meanwhile we are going to introduce you to some of the available tools. 

There are several software packages available to help your study of differ- 
ential equations. These tend to fall into two categories. On one side are 
programs which are visual and easy to use, with the focus being on using 
the computer to show you what some of the possibilities are. You point and 
click and try things, and see a lot without working too hard. On the other 
side are programs in which you have to do some active programing yourself. 
While they often have a visual component, you have to get your hands dirty 
to get anything out of these. The result is that you get a more concrete un- 
derstanding of various fundamentals. Both kinds have advantages. It seems 
that the best situation is not to choose one exclusively, but to have and use 
both kinds. 

In the easy visual category there are java applets available at 

http : / /www . math . Cornell . edu/ ~bterrell/ de 

and 

http : / /math . rice . edu/ df ield/ df pp . html 
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Probably the earliest userfriendly DE software was MacMath, by John Hub- 
bard. MacMath was the inspiration for de. 

The other approach is to do some programing in any of several available lan- 
guages. These include mat lab, its free counterpart octave, and the freeware 
program xpp. 

There are also java applets on partial differential equations. These are called 
Heat Equation ID, Heat Equation 2D, Wave Equation ID, and Wave Equa- 
tion 2D, and are available from 

http : / /www . math . Cornell . edu/ ~bterrell 
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Figure 7.1 This is the slope field for the Logistic Equation as made by 
dfield. Question: Do these curves really touch? Read the Theorem again 
if you're not sure. 

Problems 

1. Find out how to download and use some of the programs mentioned in this lecture, try 
a few simple things, and read some of the online help which they contain. 

2. Answer the question in the caption for Figure 7.1. 



24 



8 Linear First Order Equations 

Today: Applications and a solution method for first order linear equations. 
Today we consider equations like ^ + ax = b or 

x + ax = b 

These are called first order linear equations. "Linear" is here a more-or-less 
archaic use of the word and means that x and x' only occur to the first 
power. We will say later what the word "linear" means in a more modern 
sense. The coefficients a and b can be constants or functions of t. You 
should check that when a and b are constant, the solutions always look 
like x(t) = c\e~ at + C2, where c\ is an arbitrary constant but C2 has to be 
something specific. When I say "check", it doesn't mean derive it, but just 
substitute into the equation and verify. Make sure you see how similar this 
is to the case when 6 = 0. 

Practice: Solve x + 2x = 3. 

Among these equations is the very first one we ever did, y' = .028y, which 
we solved by separation of variables. Some, but not all of these equations 
can be solved the same way. The simplest one of them is x' = which is so 
easy as not to be useful. The next simplest is x' = b and this one is very 
important because it holds a key to solving all the others. We call this one 
the "easy kind" of differential equation because it does not involve x at all 
and can be solved just by integrating: 




Wouldn't it be nice if all differential equations could be solved so easily? 

Maybe some can. Somebody once noticed that all these equations can be 
put into the form of the easy kind, if you will only multiply the equation by 
a cleverly chosen "integrating factor" . For example, consider again 

x + 2x = 3 

Multiply by e 2t and you will get 
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e 2 V + 2e 2t x = 3e 2t 



Then, and this is the cool part, recognise the left hand side as a total deriva- 
tive, like you have in the easy kind: 



This uses the product rule for derivatives. Now all you have to do is integrate 
and you are done. 



What do you think the integrating factor should be for x' — 1.3x = cos(t)? 
In this case, which will be a homework problem, you still have an integral 
to do which is a little harder, and you probably have to do integration by 
parts. There is a general formula for the integrating factor which is not 
too important, but here it is. For the equation above, x' + ax = b, an 
integrating factor is e^ adt . This always works, at least within the limits of 
your ability to do the integrals. As for where this formula comes from, we 
have devised a homework problem so that you can answer that yourself. As 
to the deeper question of how anybody ever thought up such a scheme as 
integrating factors in the first place, that is much harder to answer. 

Now let's look at some applications of these equations. 
8.1 Newton's Law of Cooling 

Newton's law of cooling is the statement that the exponential growth equa- 
tion x' = kx applies sometimes to the temperature of an object, provided 
that x is taken to mean the difference in temperature between the object 
and its surroundings. 

For example suppose we have placed a 100 degree pizza in a 600 degree 
oven. We let x(t) be the pizza temperature at time t, minus 600. This 
makes x negative, while x' is certainly positive because the pizza is heating 
up. Therefore the constant k which occurs must be a negative number. 
The solution to the equation does not require an integrating factor; it is 
x(t) = Ce kt , and C = 100 — 600 as we have done previously. Therefore the 



(e 2t x)' = 3e 



,2t 




so 



x{t) = \ + Ce-- 
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pizza temperature is 600 + x(t) = 600 — 500e . We don't have any way to 
get k using the information given. It would suffice though, to be told that 
after the pizza has been in the oven for 15 minutes, its temperature is 583 
degrees. This says that 583 = 600 — 500e 15fc . So we can solve for k and 
then answer any questions about temperature at other times. Notice that 
this law of cooling for the pizza can be written as p' — kp = 600k where p{t) 
is the pizza temperature, which is a first order linear equation. A different 
environment arrives if we move the pizza to the 80 degree kitchen. A plot 
of the temperature history under such conditions is in Fig 4.1. It is not 
claimed that the solution is differentiable at t = 6, when the environment 
changes. 




Figure 8.1 The pizza heats and then cools. Note the trick used to change 
the environment from 600 degrees to 80 degrees at time 6: the equation was 
written as p'=-p+a-b*sign(t-6) , and a and b selected to achieve the 600 
and 80. (figure by df ield) 

There is another use of Newton's law of cooling which is not so pleasant to 
discuss, but it is real. The object whose temperature we are interested in 
is sometimes of interest because the police want to know the last time it 
was alive. The body temperature of a murder victim begins at 37 C, and 
decreases at a rate dependent on the environment. A 30 degree road is one 
environment, and a 10 degree morgue table is another environment. Part 
of the history is known, namely that which occurs once the body has been 
found. Measurements of time and temperature can be used to find k. Then 



27 



that part of the history from the unknown time of death to the known time 
of discovery has to be reconstructed using Newton's law of cooling. This 
involves functions not unlike those in Figure 8.1. It should be clear that this 
involves estimates and inaccuracies, but is useful, nonetheless. 

8.2 Investments 

Moving to a more uplifting subject, our bank account equation y' = .028y 
can be made more realistic and interesting. Suppose we make withdrawals 
at a rate of $3500 per year. This can be included in the equation as a 
negative influence on the rate of change. 

y = .028y - 3500 

Again we see a first order linear equation. But the equation is good for more 
than an idealised bank account. Suppose you buy a car at 2.8% financing, 
paying $3500 per year. Now loosen up your point of view and imagine 
what the bank sees. From the point of view of the bank, they just invested 
a certain amount in you, at 2.8% interest, and the balance decreases by 
"withdrawals" of about $3500 per year. 

So the same equation describes two apparently different kinds of investments. 
In fact it is probably more realistic as a model of the loan than as a bank 
account. By changing the —3500 to +3500 you can model still other kinds 
of investments. 

Example: A car is bought using the loan as described above. If the loan 
is to be paid off in 6 years, what price can we afford? 
The price is y(0). We need 

4y' = .0281/ - 3500 with y(6) = 0. 

Calculate e - 028t y'-.028e- 028t y = -3500e- 028i , (e- 02at y)' = -3500e- 028t , 
e -.028t y = (3500/.028)e- 028t +c?/ = (3500/.028) + ce' 028t = 125000 + ce' 028i , 
2/(6) = = 125000 + ce- 028(6) . This implies c = -105669, 2/(0) = 125000 - 
105669 = 19331 

Problems 

1. Solve x' — 1.3x = cost. 

2. This problem shows where the formula e$ ait for the integrating factor comes from. 
Suppose we have the idea to multiply x' + ax = b by a factor /, so that the result of the 
multiplication is (fx)' = fb. Show that for this plan to work, you will need to require 
that /' = af. Deduce that e^ adt will be a suitable choice for /. 
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3. The temperature of an apple pie is recorded as a function of time. It begins in the 
oven at 450 degrees, and is moved to an 80 degree kitchen. Later it is moved to a 40 
degree refrigerator, and finally back to the 80 degree kitchen. Make a sketch somewhat 
like Figure 8.1, which shows qualitatively the temperature history of the pie. 

4. A body is found at midnight, on a night when the air temperature is 16 degrees C. Its 
temperature is 32 degrees, and after another hour, its temperature has gone down to 30.5 
degrees. Estimate the time of death. 

5. Sara's employer contributes $3000 per year to a retirement fund, which earns 3% 
interest. Set up and solve an initial value problem to model the balance in her fund, if it 
began with $0 when she was hired. How much money will she have after 20 years? 

6. Show that the change of variables x = 1/y converts the logistic equation y = .028(y — 
y 2 ) of Lecture 4 to the first order equation x' = — .028(a; — 1), and figure out a philosophy 
for why this might hold. 

7. A rectangular tank measures 2 meters east-west by 3 meters north-south and contains 
water of depth x(t) meters, where t is measured in seconds. One pump pours water in 
at the rate of 0.05 [m 3 /sec] and a second variable pump draws water out at the rate of 
0.07 + 0.02cos(cjt) [m 3 /sec]. The variable pump has period 1 hour. Set up a differential 
equation for the water depth, including the correct value of ui. 




Figure 8.2 Here x(i) is the length of a line of people waiting to buy tickets. 
Is the rate of change proportional to the amount present? Does the ticket 
seller work twice as fast when the line is twice as long? 
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9 [oiler's] Euler's Numerical Method 



Today: A numerical method for solving differential equations either by 
hand or on the computer, several ways to run it, and how your calculator 
works. 

Today we return to one of the first questions we asked. "If your bank balance 
y(t) is $2000 now, and dy/dt = .02% so that its rate of change is $56 per 
year now, about how much will you have in one year?" Hopefully you guess 
that $2056 is a reasonable first approximation, and then realize that as soon 
as the balance grows even a little, the rate of change goes up too. The 
answer is therefore somewhat more than $2056. 

The reasoning which lead you to $2056 can be formalised as follows. We 
consider 

x' = f(x,t) 
x(to) = x 

Choose a "stepsize" h and look at the points t\ = to + h, ^2 = ^0 + 2/i, 
etc. We plan to calculate values x n which are intended to approximate the 
true values of the solution x{t n ) at those time points. The method relies on 
knowing the definition of the derivative 

x'(t) = lim X(t + /l) -" (t) 
h-^Q h 

We make the approximation 

X ' ' = h 

Then the differential equation is approximated by the difference equation 

= f i,Xn-,tn) 

Example: With h = 1 the bank account equation becomes 

Vn+X ~ Vn = .028y„ 

or 

y n +i = Vn + -028hy„ 



30 



This leads to 3/1 = yo + .028/ij/o = I.O283/0 = 2056 as we did in the first 
place. For a better approximation we may take h = .2, but then 5 steps are 
required to reach the one-year mark. We calculate successively 



yi 


= yo + .028hy 


= 1.0056y 


= 2011.200000 


1)2 


= 3/1 + .028h yi 


= 1.0056j/i 


= 2022.462720 


ys 


= y 2 + .028hy 2 


= 1.00 563/2 


= 2033.788511 




= j/ 3 + .028/ii/3 


= 1.00 563/3 


= 2045.177726 


ys 


= 3/4 + .028ftj/4 


= 1.00 563/4 


= 2056.630721 



Look, you get more money if you calculate more accurately! 

Here the bank has in effect calculated interest 5 times during the year. 
"Continuously Compounded" interest means taking h so close to that you 
are in the limiting situation of calculus. 

Now at this point we happen to know the answer to this particular problem. 
It is y(l) = 2000e' 028 ( 1 ) = 2056.79 1 369.... Continuous compounding gets 
you the most money. Usually we do not have such formulas for solutions, 
and then we have to use this or some other numerical method. 

This method is called Euler's Method, in honour of Leonard Euler, a math- 
ematician of the 18th century. He was great, really. The collected works 
of most famous scientists typically fill a few books. To see all of his in the 
library you will have to look over several shelves. This is more impressive 
still when you find out that though he was Swiss, he worked in Russia, wrote 
in Latin, and was in later life blind. Do you know what "e" stands for, as 
in 2.718...? "exponential" maybe? It stands for Euler. He also invented 
some things which go by other people's names. So show some respect, and 
pronounce his name correctly, "oiler" . 

Example: Now we'll do one for which the answer is not as easily known 
ahead of time. Since the calculations can be tedious, this is a good time to 
use the computer. We will set up the problem by hand, and solve it two 
ways on the machine. Assume that p(t) is the proportion of a population 
which carries but is not affected by a certain disease virus, initially 8%. The 
rate of change is influenced by two factors. First, each year about 5% of the 
carriers get sick, so are no longer counted in p. Second, the number of new 
carriers each year is about .02 of the population but varies a lot seasonally. 
The differential equation is 

p = -.05p+ .02(1 + sin(27rf)) 
p(0) = .08 
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Euler's approximation is to choose a stepsize h, put t n = nh, po = .08, and 
approximate p(t n ) by p n where 

Pn+i = Pn + h(-.05p„ + .02(1 + sin(t„))) 

To compute in octave or matlab, use any text editor to make a function file 
to calculate the right hand side: 



°/,file myfunc.m 

function rate = myfunc(t,p) 

rate = -.05*p + . 02* (l+sin(2*pi*t) ) ; 

"/.end file myfunc.m 



Then make another file to drive the computation: 

'/.file oiler. m 
h = .01; 
p(l) = .08; 
for n = 1:5999 

p(n+l) = p(n) + h*myfunc(n*h,p(n)) ; 
end 

plot(0:h:60-h,p) 
"/.end file oiler. m 

Then at the octave prompt, give the command oiler. You should get a graph like the 
following. Note that we used 6000 steps of size .01, so we have computed for < t < 60. 
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Figure 9.1 You can see the seasonal variation plainly, and there appears to 
be a trend to level off. This is a dangerous disease, apparently, (figure made 
by the plot commands shown in the text) 
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For additional practice you should vary h and the equation to see what happens. 

There are also more sophisticated methods than Eider's. One of them is 
built into matlab under the name ode23 and you are encouraged to use it 
even though we aren't studying the internals of how it works. You should 
type help ode23 in matlab to get information on it. To solve the carrier 
problem using ode23, you may proceed as follows. Use the same myfunc.m 
file as before. Then in Matlab give the commands 

[t,p] = ode23( , myfunc' ,0,60, .08) ; 
plot(t,p) 

Your results should be about the same as before, but generally ode23 will 
give more accuracy than Euler's method. There is also an ode45. We won't 
study error estimates in this course. However we will show one more example 
to convince you that these computations come close to things you already 
know. Look again at the simple equation x' = x, with x(0) = 1. You know 
the solution to this by now, right? Euler's method with step h gives 

^n+l — %n ~t~ hx n 

This implies that 

x\ = (1 + h)xQ = 1 + h 

X 2 = (l + h)x! = {l + hf 

x n = {l + h) n 

Thus to get an approximation for x(l) = e in n steps, we put h = 1/n and 
receive 

e = (l + -) n 
n 

Let's see if this looks right. With n = 2 we get (3/2) 2 = 9/4 = 2.25. With 
n = 6 and some arithmetic we get (7/6) 6 = 2.521626, and so forth. The point 
is not the accuracy right now, but just the fact that these calculations can be 
done without a scientific calculator. You can even go to the grocery store, 
find one of those calculators that only does H — */, and use it to compute 
important things. 

Did you ever wonder how your scientific calculator works? Sometimes people 
think all the answers are stored in there somewhere. But really it uses ideas 
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and methods like the ones here to calculate many things based only on H — */. 
Isn't that nice? 



Example: We'll estimate some cube roots by starting with a differential 
equation for x(t) = t 1/:i . Then x(l) = 1 and x'(t) = \t~ 2/3 . These give the 
differential equation 

1 



X ~ 3x2 



Then Euler's method says x n +i — x n + g^r, and we will use xo = 1, h = .1: 
xi = 1 + ^ = 1-033333... 



Therefore (1.1) 1/3 = 1.0333. 



x 2 = xi + = 1.0645 

oX\ 

Therefore (1.2) 1 / 3 = 1.0645 etc. For better accuracy, h can be decreased. 
Problems 

1. What does Euler's method give for \/2, if you approximate it by setting x(t) = y/t and 
solving 

x = — with x(l) — 1 
Use 1, 2, and 4 steps, i.e., h = 1, .5, .25 respectively. 

2. Solve 

y — (cosy) 2 with i/(0) = 
for < t < 3 by modifying the oiler. m and myfunc.m files given in the text. 

3. Solve the differential equation in problem 2 by separation of variables. 

4. Redo problem 2 using ode45 and your modified myfunc.m, as in the text. 

5. Run df ield on the equation of problem 2. Estimate the value of j/(3) from the graph. 

6. Compare your answers to problems 2 and 3. Is it true that you just computed tan _1 (3) 
using only H — */ and cosine? Figure out a way to compute tan _1 (3) using only H — */. 

7. Solve the carrier equation p' = — .05p + .02(1 + sin(27rf)) using the integrating factor 
method of Lecture 8. The integral is pretty hard, but you can do it. Predict from your 
solution, the proportion of the population at which the number of carriers "levels off" 
after a long time, remembering from Figure 8.1 that there will probably continue to be 
fluctuations about this value. Does your number seem to agree with the picture? Note 
that your number is probably more accurate than the computer's. 

8. Newton's law of cooling looks like u' — —au when the surroundings are at temperature 
0. This is sometimes replaced by the Stephan-Boltzmann law v! = —bu 4 , if the heat is 
radiated away rather than conducted away. Suppose the constants a, b are adjusted so 
that the two rates are the same at some temperature, say 10 Kelvin. Which of these laws 
predicts faster cooling when u < 10? u > 10? 
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10 Second Order Linear Equations 



Today: 2nd order equations and some applications. The characteristic 
equation. Conservation laws. 

The prototype for today's subject is x" = —x. You know the solutions to 
this already, though you may not realize it. Think about the functions and 
derivatives you know from calculus. In fact, here is a good method for any 
differential equation, not just this one. Make a list of the functions you 
know, starting with the very simplest. Your list might be 



1 

c 

t 

t n 
e* 
cos(i) 

Now run down the list trying things in the differential equation. In x" = —x 
try 0. Well! what do you know? It works. The next few don't work. Then 
cos(i) works. Also sin(t) works. Frequently, as here, you don't need to use a 
very long list before finding something. As it happens, cos(t) and sin(t) are 
not the only solutions to x" = —x. You wouldn't think of it right away, but 
2cos(t) — 5sin(£) also works, and in fact any c\ cos(t) + C2 sin(t) is a solution. 

Practice: Find similar solutions to x" = —9a:. 

The equations x" = cx happen to occur frequently enough that you should 
know all their solutions. Try to guess all solutions of x" = and x" = 25a; 
for practice. 

10.1 Linearity 

We will consider equations of the forms 

ax" + bx + cx = 
ax" + bx + cx = dcos(ft) 
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These are called second order linear equations, and the first one is called 
homogeneous. The strict meaning of the word 'linear' is being stretched 
somewhat here, but the terminology is traditional. Really an operation L 
which may be applied to functions, such as L(x) = ax" + bx' + cx, is linear 
if it is true that L(x + y) = L(x) + L(y) and L(kx) = kL(x). This does hold 
here, and the consequence is that if x and y are both solutions to L{x) = 0, 
then linear combinations like 2x + by are also solutions. That is true for 
the homogeneous equation because we have L(2x + by) = L{2x) + L(5y) = 
2L(x) + 5L(y) = 0. However, the nonhomogeneous second equation does 
not have this property. You can see this yourself. Suppose x and y are 
solutions to the same equation, say x" + 3x = 2 and y" + 3y = 2. Then 
we want to see whether the sum s = x + y is also a solution. We have 
s" + 3s = x" + y" + 3x + 3y = 2 + 2 = 4. This is a different equation, 
since we got s" + 3s = 4, not s" + 3s = 2. So it is a bit of a stretch to call 
the nonhomogeneous equation 'linear', but that is the standard terminology. 
You should also be sure you understand that for an equation like x" ' +x 2 = 0, 
even with on the right side, the sum of solutions is not usually a solution. 

10.2 The Characteristic Equation 

Next we look at a method for solving the homogeneous second order linear 
equations. Some of the nonhomogeneous ones will be done in Lecture 13. 
The motivation for this method is the observation that exponential functions 
have appeared several times in the equations which we have been able to 
solve. Let us try looking down our list of functions until we come to the 
ones like e rt . Trying x = e rt in 

ax" + bx' + cx = 

we find ar 2 e rt + bre rt + ce rt = (ar 2 + br + c)e rt This will be zero only if 

ar 2 + br + c = 

since the exponential is never 0. This is called the characteristic equation. 
For example, the characteristic equation of x" + Ax 1 — 3x = is r 2 + 4r — 
3 = 0. See the similarity? We converted a differential equation to an 
algebraic equation that looks abstractly similar. Since we ended up with 
quadratic equations, we will get two roots usually, hence two solutions to 
our differential equation. 
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Example: x" + Ax' — 3x = has characteristic equation r 2 + 4r — 3 = 0. 
Using the quadratic formula, the roots are r = — 2 ± \fj. So we have found 
two solutions x = e^ 2- ^'* and x = e l '~ 2+ ^ t . By linearity we can form 
linear combinations 

Example: You already know x" = —x very well, right? But our method 
gives the characteristic equation r 2 + l = 0. This does not have real solutions. 
In Lecture 12 we'll see how the use of complex numbers at this stage will 
unify the sines and cosines we know, with the exponential functions which 
our characteristic method produces. 

Example: The weirdest thing which can happen is the case in which the 
quadratic equation has a repeated root. For example, (r + 5) 2 — r 2 + lOr + 
25 = corresponds to x" + 10x' + 25a; = 0. The only root is —5, so we 
only get one solution x(t) — e~ 5i , and multiples of that. If these were the 
only solutions, we could solve for say, the initial conditions x(0) — 3 and 
x'(0) = —15, but could not solve for x(0) = 3 and x'(0) = —14. As it turns 
out, we have not made a general enough guess, and have to go further down 
our list of functions, so to speak, until we get to things like te rt . These work. 
So we end up with solutions x(t) = (ci + C2t)e~ 5t for this problem. This is a 
special case, and not too important, but it takes more work than the others 
for some reason. 



Example: 

x" = —5a; 
x(0) = 1 
x'(0) = 2 

General solution is x(t) — c\ cos(\/5£) + C2 sin(\/5t). Use initial conditions 
x(0) = ci = 1, x'(0) = V5c 2 = 2,c 2 = ^. So 

2 

x(t) — cos(v / 5t) H — -j= sin(\/5i) 



Example: 

y" = 4y 
2/(0) = 1 
y'(0) = 2 

Characteristic equation is r 2 = 4, r = ±2, y — c\e 2t + C2e~ 2t , y(0) = 
C-L + C2 = 1 and y'(0) = 2cj - 2c 2 = 2. Then eliminate c 2 : 2y(0) + y'(0) = 
4ci = 4, ci = 1, C2 = 1 — ci = 0, y(t) — e 2i , you check it. 
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10.3 Conservation laws, part 2 



Sometimes a second order equation can be integrated once to yield a first 
order equation. Good things happen when this is possible, since we have 
many tools to use on first order equations, and because usually the resulting 
equation can be interpreted as conservation of energy or some similar thing. 

For example, let's pretend that we don't know how to solve the equation 
x" = —x. You can try to integrate this equation with respect to t. Look 
what happens: 

/***=-/** 

You can do the left side, getting x', but what happens on the right? You 
can't do it can you? The right side doesn't integrate to xt or something 
like that does it? I'm asking all these questions because this is a common 
place to make mistakes. Really, you just can't do the integral on the right, 
because the integrand x(t) is not yet known. But there is a way to fix this 
problem. It is very clever because it seems at first to complicate the problem 
rather than simplifying it. Multiply the equation x" = —x by x' and see 
what happens. You get 

x'x" = —xx 1 

Now stare at that for a while, and see if you can integrate it now. 

Jx'x"dt = -jxx' d t 

You may notice, if you look at this long enough, that according to the chain 
rule you can now integrate both sides. That is because xx' is the derivative 
of ^(x) 2 , and x'x" is the derivative of ^(x') 2 . So integrating, you get 

\{x'? = c-\{xf 

Now there are several things to observe about this. First we don't have a 
second order equation any more, but a first order instead. Second we can 
make a picture for this equation in the (x, x') plane, and it just says that 
solutions to our problem go around circles in this plane, since we just got the 
equation of a circle whose size depends on the constant c. Third, there is a 
physical interpretation for the first order equation, which is conservation of 
energy. Not everyone who takes this course has had physics yet, so the next 
bit of this discussion is optional, and may not make too much sense if you 
have not seen an introduction to mechanics. Conservation of energy means 
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the following, x is the position and x' the velocity of a vibrating object. 
Peek ahead to Figure 13.1 if you like, to see what it is. The energy of this 
object is in two forms, which are kinetic energy and potential energy. The 
kinetic energy "^mu 2 " is the (x') 2 term. The potential energy u ^kx 2v is the 
x 2 term. So what is c? It is the total energy of the oscillator, and the going 
around in circles is physically the fact that when something swings back and 
forth, the energy is periodically transferred from to potential to kinetic and 
back. The child on a swing has a lot of kinetic energy at the bottom of the 
path moving fast, and all that becomes potential energy by the top of the 
swing when he is instantaneously not moving, and the cycle repeats. 

This kind of calculation, finding the integral of the equation or some multiple 
of the equation, is called finding a conservation law. It cannot always be 
used, but is a powerful thing to know about. 

Here is an example of the power. It is a little more advanced than the rest 
of the lecture, but I bet you will like it. 

Theorem There is no other realvalued solution to x" = —x than the ones 
you already know about. 

You probably wondered whether anything besides the sine and cosine had 
that property. Of course there are the linear combinations of those. But 
maybe we just aren't smart enough to figure out others. The Theorem says 
no: that's all there are. 

PROOF Suppose the initial values are x(0) = a and x'(0) = b, and we 
write down the answer x(t) that we know how to do. Then suppose your 
friend claims there is a second answer to the problem, called y(t). Set 
u(t) — x(t) — y(t) for the difference which we hope to prove is 0. Then 
u" — —u. We know from the conservation law idea that then 

What is c? The initial values of u are 0, so c = 0. Do you see why that 
makes u identically 0? 

Problems 

1. Suppose x'{ + xi 2 = and x'{ + X2 2 = 0, and set s = xi + xi- Calculate s" + s 2 to see 
why the sum of solutions is not usually a solution. 

2. Suppose xi and X2 both solve 

x + x = 
x'(0) = 2 

and set s = x\ + x 2 . Find out why s is not a solution. 
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3. Suppose xi and £2 both solve 



x" + x = cos 8t 
x'(0) = 



and set s = x\ + xi. Is s a solution? 



4. Solve 



W - 3y' -y 
2/(0) 

y'(o) 



o 







1 



5. Solve the same equation as in problem 4, but 



with initial conditions 



2/(0) = 
2/(0) = 







1 



6. Solve the same equation as in problem 4, but 



with initial conditions 



2/(0) = 
2/'(0) = 



2 



3 



7. Can you see that the answer to problem 6 should be 2 times the answer to problem 4 
plus 3 times the answer to problem 5? 

8. Find a conservation law for the equation x" + x [i = 0. 

9. Do you think there are any conservation laws for x" + x' + x = 0? 

10. What's rong with this? y" + y 2 = 0, r 2 + 1 = 0, r = -1, y = e~* + c. There are at 
least 3 errors. 
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11 The Heat Equation 



Today: We now have enough information to make a good start on the 
partial differential equation of heat conduction. If you want to save PDE for 
later, you could read this lecture after Lecture 28. 

The heat equation looks like this: 

du ,9 2 u 
dt = dx* 

usually abbreviated ut = ku XXl and it describes the temperature of an object 
in which heat energy is allowed to flow by conduction in the x direction. 
Heat can also move with a flowing fluid but that, known as convection, is 
not what this is about. Heat can also move as radiation, like the warming 
you feel near a fire even though the air itself may be cold. It is not about 
that either. It can describe the changing temperature along a metal bar or 
through a thick wall. The k is a material constant. For now, take k = 1. 
Later we give a derivation of this from physical principles. Right now, lets 
see how close we can come to solving it. 

Try the separation idea. Look for solutions of the form 

u(x,t) = X(x)T(t) 

This does not mean we think they are all like this. We hope that some are. 
Substituting, we need 

T'(t)X(x) = T(t)X"(x) 

Now divide by XT, assuming for now that this is not 0. 

Of course, later we might find that sometimes it is 0, but by then we will have 
made an independent check, without dividing anything, that our answer 
really is a correct answer. So right now we are just exploring, not really 
proving anything. We get 

~T ~ ~X~ 

Again we see that since the left side only depends on t, and the right does 
not depend on t, that neither side depends on t, or on x. They are the same 
constant, call it — c. We must then solve the ordinary differential equations 

T' = -cT and X" = -cX 
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Observe, two good things happened simultaneously. First the PDE broke 
into two ODEs. Second the ODEs turned out to be some that we already 
know how to solve. 



Example: Suppose c = 1. We know how to solve T" = -T, and X" = -X. 
We know solutions to those. T(t) — Cie - *, and X(x) = C2 cos(:r) + C3 sin(a;). 
Thus some solutions to the heat equation ought to be 



u(x, t) 



e cos(a;) 



u(x,t) — 3e sin(a;) 
u(x,t) — e~ (6cos(a;) — 5.3sin(a:)) 



Now, we must check those. For one thing, they are in a few places and 
maybe that matters. Or maybe we overlooked something else in our sepa- 
ration argument. 

Practice: Plug those into the heat equation and verify that they really do 
work. 

Problems 

1. Suppose c = 2 instead of 1 in our separated ODEs 



What solutions T, X, u do you find? 

2. What if c is some other positive number besides 1 or 2? What solutions T, X, u do 
you find? 

3. What happens if c = 0? What solutions T, X, u do you find? 

4. Sketch a graph of the temperature for a few times, in the solution u(x, t) — 3e~* sin(x) 



5. What happens if c is a negative number, say — 1? Sketch a graph of the temperature 
for a few times. Could that u be the temperature of an object? There is something wrong 
with that as a object temperature in most cases, isn't there? So we don't use negative c 
usually. 

6. Find some solutions to the heat equation 



T' = - 



cT 



and X" — —cX 



above. 



ut = ku. 



'XX 



where the material constant k has been put back in. 
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12 Complex Numbers and the Characteristic Equa- 
tion 



Today: The complex numbers, and a motivation for the definition of e a+ '. 
More second order equations. 

Complex numbers are expressions of the form a + bi where a and b are real 
numbers. You add, subtract, and multiply them just the way you think you 
do, except that i 2 = —1. So for example, 

(2.5 + 3i) 2 = 6.25 + 2(2.5)(3i) + 9i 2 = -3.25 + 15* 

If you plot these points on a plane, plotting the point (x, y) for the complex 
number x + yi, you will see that the angle from the positive x axis to 2.5 + 3i 
gets doubled when you square, and the length gets squared. Addition and 
multiplication are in fact both very geometric, as you can see from the 
figures. 



5 



3 








Figure 12.1 Addition of complex numbers is the same as for vectors. Mul- 
tiplication adds the angles while multiplying the lengths. The left picture 
illustrates the sum of 3 + i and 1 + 2i, while the right is for the product, 
(figure made by the compass command in matlab) 

Practice: Use the geometric interpretation of multiplication to figure out 
a square root of i. 
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Division of complex numbers is best accomplished by using this formula for 
reciprocals: 

1 a — bi 

a + bi a 2 + b 2 



Practice: Verify that this reciprocal formula is correct, i.e. that using it 

You are now equiped to do arithmetic, such as solving 2 — (2 — 5i)z = 3/i for 
z. To do algebra correctly we should point out that two complex numbers 
are equal by definition when the real and imaginary parts are equal. If 
z = a + bi and a and b are real, then we write re(z) = a and im(z) = b. (not 
bi.) Warning: If you know that a + bi = c + di you cannot conclude that 
a = c and b = d. For example, 1 + Oi = + {—i)i- However, if you know in 
addition that a, b, c, and d are real, then you can conclude that a = c and 
b = d. 

The quadratic formula for solving ax 2 + bx + c = is still good if everything 
in sight is complex, but taking square roots can be tedious. That is why we 
practiced it above. 

There is another important construction which is motivated by algebra and 
differential equations, and very useful. That is raising e to a complex power. 
We set e^ s+t '^ = e s e tl by analogy with known properties of real exponentials, 
but this still requires a definition of e tl . We claim that the only reasonable 
choice is cos(t) + isin(i). The reason is as follows. The whole process 
of solving second order equations by the characteristic equation method 
depends on the formula 



d{e rt ) 



re 



dt 

Let's require that this hold also when r = i. Writing e lt = f(t) + ig(t) this 
requires 

/' + ig' = i(f + ig) =-g + if 

so /' = —g and g' = f. These should be solved using the initial conditions 
e° = 1 = /(0) + ig(0). These give /(0) = 1, /'(0) = and /" = -/. The 
only solution is f(t) = cos(t) and g(t) = sin(t). Therefore our definition 
becomes 

e (s+ti) = e s( cos ( t ) + j s i n (t)) 

Now it turns out that we actually get ( ,. ' = re for all complex r. You 
should verify this for practice, using our definition. This is a case of pulling 
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oneself up by one's own bootstraps: we solve a simple differential equation 
to make a definition of e to a complex power, so that we can use that to 
solve more complicated differential equations. 

Now we can do more examples. 

Example: x" + 2x' + 3x = The characteristic equation, which works 
by all that has been done above, is r 2 + 2r + 3 = 0. The solutions are 
r = — 1 ± y/—2 = — 1 ± iV2. So there are solutions to the differential 
equation of the form 

x(t) = c 1 e ( - 1+i ^ )t + c 2 e ( - 1 - < ^ )t 

Given some initial conditions, we can find the constants. For instance if we 
want to have x(0) = 5 and x'(0) = then we must solve 

x(0) —5 — Cl + C2 

x'(0) = = (-l + iv / 2)ci + (-l-iv / 2)c 2 
The sum of these gives 5 = iv2(ci —02). So 

5 5i 

Cl = C2 H t= — C2 t= 

The first equation then gives 

5 = 2c 2 -^ 

Finally c 2 = (5/2)(l - and ci = 5 - c 2 = (5/2)(-l - ^). So 
5 / . 1 \ (-l+iv^t , 5 / J_\ ( _i_iV2)t 



You may notice that the answer, besides being messy, looks very complex. 
Now the last time you went into the lab, all the instruments were probably 
reading out real numbers, weren't they? So it is common to rewrite our 
solution in a way which shows that it really is real after all. To do this, the 
following new idea is required: If x(t) is a complex solution to a second order 
linear differential equation with real coefficients then the real and imaginary 
parts of x are also solutions! For example, e 3tt is a solution to x" = —9x. 
The real and imaginary parts are respectively cos(3t) and sin(3t), and these 
are certainly solutions also. To see why this works in general, suppose that 
x = u + iv solves ax" + bx' + cx = 0. This says that a(u" + iv") + b(u' + iv') + 
c(u + iv) = (au" + bu' + cu) + i(av" + bv' + cv) = Assuming that a, b, c, u, 
and v are all real, you can conclude that au" + bu' + cu and av" + bv' + cv are 
also 0. But note that if some of the coefficients had not been real, or if the 
equation had been something like x" + x 2 = then the argument would not 
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have worked. Equations with explicit complex coefficients sometimes occur 
in physics, in spite of what we all noticed above about lab instruments. 
Schrodinger's Equation is an example. In such cases you can't simply use 
the real and imaginary parts. 

So let's rework the previous example slightly. After finding the values of 
r above, we wrote down among other things a complex solution x(t) = 
e (-iW2)t = e -*(cos( v / 2t) +ism(V2t)) Using the real and imaginary parts, 
we form solutions 

x{t) = aie^ 1 cos(v / 2t) + a 2 e~* sin(v / 2t) 

If we apply the same initial conditions again, we need x(0) = 5 = aq and 
x'(0) = = — cti + \/2ci2- Solving these we end with 

x(t) = 5e~* cos(v / 2t) - -^e _< sin( v / 2t) 

It is important to realise that this is equal to the complex-looking solution 
we found previously. 



Example: Solve for a: 

?-(2 + i)o = 7 



"(2+0 



7+3i _ _/ 7 , o -\ 2-i _ -14-K? + 7i+:ii^ _ -17+i 



Example: Solve 

y" + 2y + 2y = 

The chara cter istic equation is r + 2r + 2 = 0, so by the quadratic formula 
r = - 2± f^ = -l±i. One solution is e ( - 1+<)< = e _t (cos(t) + isin(t)). 
Taking the real and imaginary parts, the solution is 

y(t) — cie - ' cos(£) + c 2 e~ J sin(t) 

Don't forget to check it. 



Example: 

x" + 2a;' + a; = 

r 2 + 2r + 1 = (r + l) 2 = 0, r = -1,-1, a; = (ci + c 2 t)e _t , Check: a;' = 
c 2 e~* - (ci + c 2 t)e^ t = c 2 e~* - x. x" = -c 2 e~* - x' — -c 2 e~' - (c 2 e~' - x) 
= -2c 2 e~' + x, x" + 2x + x = -2c 2 e~' + x + 2(c 2 e~* - x) + x = 
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Example: 

x" + 2x — x — 

r 2 + 2r - 1 = 0, r = -1 ± VTTl, x(t) = Cl e<- 1 -^*+ ae^ 1 ^* 
Problems 

1. Solve r 2 - 6r + 10 = 0, r A - 6r 2 + lOr = 0, and r 2 - 6r - 10 = 0. 

2. Solve y" + 3y' + 4y = 0. 

3. Find a, if a + i = 0; if a + i = i. 

4. Sketch the graphs of the functions e _< cos(i), cos(3i), and e~ 4 * sin(3t). 

5. Verify that = re rt for r = a + bi, using the definition e (a+M)t = e at (cos(fe) + 
i sin(6t)). 

6. For each t, e' 1+l '* is a complex number, which is a point in the plane. So as t varies, a 
curve is traced in the plane. Sketch it. 

7. A polynomial r 2 + br + c = has roots r = —2 ± i. Find 6 and c. 

8. A characteristic equation r 2 + br + c — has roots r = — 2±i. What was the differential 
equation? 

9. Miai rony wtfi i/iis? x" + 4a; = 0, r 2 + 4 = 0, r = ±2i, x = d{e 2it + e" 2li ). 

12.1 The Fundamental Theorem of Algebra 

Talk about lifting yourself by your own bootstraps. A new number writ- 
ten i was invented to solve the equation x 2 + 1 = 0. Of course you have 
noticed that the people who invented it were not very happy about it: con- 
trast "complex" and "imaginary" with "real" and "rational" . But soon the 
following theorem was discovered. 

By introducing the new number i, not only can you solve x 2 + 1 = 0, but 
you can also solve at least in principle all these: x 2 + 2 = 0, x 3 — i = 0, 
x 6 - (3 - 2i)x A - x 3 + iirx - 39.1778i + 43.2 = 0, etc, etc. 

Fundamental Theorem of Algebra Let 

p(x) = a + a\x + a 2 x 2 H h a n _ix n-1 + x n 

be a polynomial of degree n > with any complex coefficients 
<2fc . Then there are complex numbers r±, ■ ■ ■ r n which are 
roots of p and p factors as 

p(x) = (x - ri)(x - r 2 ) • • • (x - r n ) 
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As a practical matter, we can handle the cases ao + x and ao + aix + x 2 easily. We all 
know the quadratic formula. There is a cubic formula that most people don't know for 
solving the n = 3 case, and a quartic formula that hardly anybody knows for the n — 4 
case, which takes about a page to write. Then something interesting happens. 

Abel's Theorem Let n > 5. Suppose you figure out every formula that 
could ever be written in terms of the coefficients au and the operations of 
addition, subtraction, multiplication, division, and extraction of roots. Then 
none of those formulas give you the roots of the polynomial. 

So: the roots exist, but if you need them then you will usually have to approximate them 
numerically. It is tough. 

Problems 

1. In the Fundamental Theorem we took the coefficient of x n to be 1 just for convenience. 
But if you don't do that the factorization must be written differently. For example, 

-6 + x + x 2 = (x-2)(x + 3) 

is correct. Figure out the number c in the case 

-30 + 5x + 5x 2 = c(x - 2)(x + 3) 

2. How must the factorization be written in general if you don't assume the coefficient of 
x n is 1? 

3. You know that in the factorization — Q + x + x' 2 = (x — 2) (a; + 3) you have —6 = (—2) (3). 
In general what is the product of the roots r\ri ■ • ■ r n in terms of parameters appearing 
in the Fundamental Theorem? 

4. You know that in the factorization —6 + x + x 2 — (x — 2)(x + 3), the x term comes 
from x — —2x + 3x when you multiply the right hand side. In general what is the sum of 
the roots n + ■ ■ ■ + r n in terms of parameters appearing in the Fundamental Theorem? 
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13 Forced Second Order 



Today: Forced equations, and how to solve them using the undetermined 
coefficients method. Natural frequency. 

Consider the equations 

x" + 2x = 
y" + 2y = 9 sin 3t 

The first one is called the homogeneous form of the second one, or the 
second is called a forced form of the first, and there are variations on this 
terminology. Mechanically what they mean is as follows. Since we know the 
solutions to the first one (don't we?) are 

x(t) = c\ cos \plt + C2 sin \p2t 

it seems reasonable that this first equation is about something vibrating 
or oscillating. It can be interpreted as a case of Newton's F = ma law, 
if you write it as —2x = lx". Here x is the position of a unit mass, x" 
is its acceleration, and there is a force — 2x which seems to oppose the 
displacement x. We call this a "spring-mass" system. It can be drawn as 
in Figure 13.1, where x is measured up. 




Figure 13.1 Unforced and forced spring-mass systems. 

The — 2x is interpreted as a spring force because it is in the direction opposite 
x: if you pull the spring 1.54 units up, then x = 1.54 and the force is —3.08, 
or 3.08 units downward. This is also a system without friction as we see 
from the fact that there are no other forces except for the spring force, and 
that the oscillation continues undiminished forever. Note that the "natural 
frequency" of this system is ^ cycles/second, in the sense that the period 
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The so-called forced equation involves an additional force, as you can see if 
you write it as — 2y — 9sin3t = ly". The picture in this case is like the right 
side of Figure 13.1. 

Now we turn to solution methods for the forced equation. We are guided 
by the physics. That is, ask yourself what will happen if you start with a 
system which wants to vibrate at a frequency of ^ , and somebody reaches 
in and shakes it at a frequency of What could happen? If you think 
about it, it is plausible that the motion might be rather complicated, but 
that a part of the motion could be at each of these frequencies. 

Let's try that. Assume 

y (t) = x(t) + Asin(3t) 

where x is the solution given above for the unforced equation. Then 

y" + 2y = x" - 3 2 Asin(3t) + 2(x(t) + Asin(3t)) = -8Asin(3t) 

We want this to equal 9sin(3i). Notice how the terms involving x dropped 
out. We need — 8A = 9 so A = —9/8. Our solution becomes 

y(t) = a cos(\/2t) + c 2 sin(\/2i) - | sin(3t) 

8 

In particular, our physics-style thinking did lead to a solution. Some people 
who take this course have had physics, and some have not had physics. In 
any case we hope you can see that a little physics, or banking, or later bio- 
logical, knowledge goes a long way in helping to set up and solve differential 
equations. 

Notice also the coefficients of the terms in this solution. There are terms 
at the natural frequency and a term at the forced frequency Jp. The 
coefficient of the forced-frequency term is fixed at — |, while the natural- 
frequency terms have arbitrary coefficients. If you have initial conditions to 
apply, now is the time. 

Example: To solve 

y" + 2y 
2/(0) 
2/(0) 

we use the solution found above 
y(t) = ci cos(v2i) 



= 9sin(3t) 
= 4 
= 5 

r- 9 
+ C2 sin(V2t) — - sin(3t) 
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We get j/(0) = 4 = ci, y'{0) = 5 = C2\/2 — 3|. The answer becomes 

5 + — 9 

y(t) = 4cos(V2t) H sin(\/2t) - - sin(3t) 

V2 8 

This is graphed in the figure 13.2. 
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Figure 13.2 Can you see two frequencies at work here? They are not the 
obvious ones! Do problems 1-3 to see what they are. 

There is a terminology that traditionally goes with this solution. One calls 
— | sin(3i) a "particular" solution to y" + 2y = 9sin(3t). It is also true that 
29.3 cos (y/2t) - | sin(3t) is a particular solution, as are many other things. 
One says that the solution is the sum of a particular solution and a solution 
of the homogeneous equation. 

Other equations can be solved by the same method. For example, it is 
possible to solve y" + y' + y = 3 + e* by guessing a particular solution of 
the form A + Be 1 . This will work. However this type of forcing is really not 
very important for us. Sines and cosines are the important things because 
they are the most practical; if you are sitting in class and happen to notice 
a fan vibration in the ceiling of the room, that is an example which can 
be modeled by one of these equations with sinusoidal forcing. The spring 
and mass correspond to some of the ceiling materials, and the sine or cosine 
forcing term represents the fan motor running at a constant speed. 

Example: 

y" + 3y' + 2y = sin(t) 

Try y = A sin(t) + B cos(t), since sine alone will not work, considering the 3y' 
term. Then y" + 3y' + 2y = {-A + 2A- 35) sin(f) + (-B + 2B + 3A) cos(i). 
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We need 

A — 3B = 1 
B + 3A = 

Substituting the second of these into the first gives A — 3(— 3^4) = f OA — 1. 
Then A = .1 and B = — 3A — —.3. Our particular solution is .lsin(t) — 
.3cos(i). For the homogeneous equation the characteristic equation is 

r 2 + 3r + 2 = (r + l)(r + 2) = 

so r — — f or r = —2. The answer is then 

y(t) = cie~* + C2e~ 2 * + .1 sin(t) — .3cos(£) 

This is shown in figure f3.3, for the case ci = 1, C2 = — 1. 

0.4 
0.3 
0.2 
0.1 


-0.1 
-0.2 
-0.3 
-0.4 
-0.5 
-0.6, 

Figure f3.3 A solution to y" + 3y' + 2y — suit. Only the forced solution 
remains after some time; the e~* and e~ 2t terms are "transient" . 

There is one special case which requires more work to solve than the others. 
This is the case which occurs when the forcing is at exactly the natural 
frequency, and is called "resonance" . Like: 

x" + Ax = cos(2t) 

Think about what happens then. If your little brother is on the swing going 
back and forth every 3 seconds, and you "help" by pushing once every 3 
seconds, pretty soon he will be going higher than he wants to go! Mathe- 
matically, if you try a particular solution ^cos^t) it won't work because 
you get when you substitute it into the left side of the equation. 

Taking these things together, it is reasonable to try something like At cos(2£) + 
Bts\n{2t). These actually work, giving solutions which grow with time. 




5 10 15 20 



52 



Practice: Verify that the solution to x"+-ix = cos(2t) is x(t) = \t sin(2t)+ 
Ci cos(2£) + c 2 sin(2£) 



Problems 



1. Use one of the addition formulas sin(a + fe) = sin(a) cos(b)+ cos(a) sin(fe) or cos(a + b) = 
cos(a) cos(6)— sin(a) sin(b) to show that any function of the form Acos(ft) + B sin(ft) 
can be written in the form C sin(/t + <f>), for some numbers C and <f>. Note that here we 
were combining sinusoids having the same frequency f/2n. 

2. Similarly to problem 1, show how a combination of two sinusoids of different frequencies 
can be written as a sum of products of sinusoids: derive for example the identity sin(8t) + 
2.3sin(7t) = 3.3sin(7.5t) cos(.5t) - 1.3cos(7.5i) sin(.5t) 

3. Combine the ideas of problems 1 and 2 to explain the two obvious frequencies in Figure 
13.2, in particular that they are not ^| and 

4. Solve x" + x + \x = 7sin6t. 

5. Solve x" + 4x — sin3t, and x" + 4x — sin2i, and discuss the difference. 

6. Show that the addition formulas for the sine and cosine are equivalent to the complex 
exponential rule e a+b = e a e b . 

7. A challenge: After a lecture involving 'frequency' one day, I noticed a crane at the 
physics building, lifting a heavy load. The main load didn't swing much, but an interme- 
diate pulley seemed to oscillate too fast for safety on such large equipment. 



Figure 13.4 A model for a crane holding a heavy load, with an intermediate 
pulley in the cable. Doesn't the main load rise twice as much as the pulley? 

Verify that the pulley rises I — I cos(#) from its lowest position. If you can work out the 
physics, try to see that the energy ignoring friction is 



Form a time derivative of that. There is an approximation sin(0) = 9 for small angles, 
which is explained better in problem 1 on page 144. Using that, derive a second order 
linear equation 8" + + l) 6 = 0. Deduce that the period of oscillation of the pulley 



What that means: you can estimate how much mass M they are lifting by just timing the 
swing of the pulley, and making a reasonable estimate of the pulley mass! 

8. What's rong with this? x" + Ax = 0, x = cos(2t) + sin(2i) + C 





is 



\l (^+l) 9 
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14 Systems of ODE's, part 1 



Today: Why would anybody want to study more than one equation at a 
time? Find out. Phase plane. 



14.1 A Chemical Engineering problem 

Jane's Candy Factory contains many things of importance to chemical en- 
gineers. One of the processing lines contains two tanks of sugar solution. 
Pure water and sugar continually enter the first one where they are mixed, 
and an equal volume of the solution flows to the second, where more sugar 
is added. The solution is taken out at the same rate from the second tank. 
The flow rates are as shown in the table. 





tank 1 


tank 2 


sugar input 


10 lb/hr 


6 lb/hr 


water input 


25 gal/hr 





solution out 


25 gal/hr 


25 gal/hr 


tank contents 


100 gal 


125 gal 


weight of sugar 


x lb 


y lb 



Jane's control system tracks the weight of sugar present in the tanks in the 
variables x and y. We can figure out the rates of change easily. The basic 
principle here is conservation of mass. 

x lib] 

x' [lb/hr] = the rate in- the rate out = 10 [lb/hr] L r J (25 [gal/hr]) 

100 [gal] 

See, you just keep an eye on the units, and everything works out. Similarly 

y' = 6 + — (25) -—(25) 
y 100 v ; 125 v ; 

Jane would like to predict the length of time required to start up this system, 
i.e., beginning with pure water in brand new tanks, how long before a "steady 
state" condition is reached (if ever)? This is an initial value problem: 

x = 10 - .25x 

y' = 6 + .25x - .2y 

x(0) = 0,y(0) = 
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Let's solve it the quickest way we can, and then we'll do more of an overview 
of systems in general. Look-the x equation is first order linear. Jane's Candy 
Factory is in luck! We know the solution is of the form x = ae~' 25 * + b. 
Using the equation and the initial condition yields x(t) = 40 — 40e - ' 25 '. 
Then the y equation reads y' = 6 + .25(40 — 40e - - 25 ') — fly which is also first 
order linear. The integrating factor is e~' 2t , which brings the y equation to 
(e- 2t y)' = 16e- 2 ' - 10e"- 05t . Then e- 2t y = 80e- 2t - 200e-- 05 * + c. The initial 
condition y(0) = gives finally c = 120 and our solution is 

x(t) = 40-40e~- 25t 

y {t) = 80 - 200e - - 25 ' + 120e"- 2 * 

The figure shows a plot of these solutions. It was made using the matlab 
commands 

t = 0:.2:30; this was a guess, that 30 hours would be enough 
x = 40 - 40*exp(-.25*t) ; 

y = 80 - 200*exp(-.25*t) + 120*exp(- . 2*t) ; 
plot(t,x); hold; plot(t,y) 




5 10 15 20 25 30 



Figure 14.1 For Jane's Candy Factory. See whether you can determine 
which curve is for the first tank, (figure made by plot commands as given 
in the text) 

Finally we can answer Jane's question about how long the plant is going to 
take to start up. As t — > oo, x(t) approaches 40 and y(t) approaches 80. 
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Using the graph, we will be within 10% of these values after about 8 hours, 
and well over half there in only 3. That kind of information could be useful, 
if you run a candy factory, or a chemical plant. 



14.2 Phase Plane 

There is an easier way to produce a plot for Jane's system, although reading 
it requires a bit of practice. The pplane utility makes the following picture. 




Figure 14.2 Jane's Candy Factory, t hidden. The limit point (40,80) corre- 
sponds to the limiting values in Figure 14.1. (figure made by pplane) 

In order to interpret this figure, one has to realise that the time is not 
shown, so that the relation between x and y may be displayed. This is 
called the "phase plane" for the system. Here we explain the idea of a phase 
plane. The vectors shown are made by using the expressions for x' and y' 
as components. This can be done for any system of the form 

x = f(x,y) 
y' = g{x,y) 

It means that we are thinking of the functions / and g as components of a 
vector field. We also think of (x(t),y(t)) as parametric equations for a curve 
in the (x, y) plane. Then from this point of view, look at what the differential 
equations say. The right-hand sides are the vector field, while the left-hand 
sides give the velocity. Consequently, the differential equation system just 
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says that the velocity of a solution curve must agree with the given vector 
field. This is the entire content of the diffential equation. You may observe 
in Figure 14.2 that the solution curves are tangent to the vectors. You do 
observe that, right? 

Practice: Choose a point in the (a;, y) plane, and work out the right hand 
sides of Jane's system for these coordinates. Then draw a vector with these 
components starting from your point (x, y). See if it is the same as the one 
pplane drew. 

Converting a System to a Higher order equation 

It is occasionally of interest to convert a system of differential equations 
into a single higher order equation. This cannot always be done. It can 
be carried out for Jane's system as follows. The y equation gives .25a; = 
y' — 6 + .2y, or x = 4y' — 24 + .8y. Substituting this into the x equation gives 
x' = Ay" - + .8y' = 10 - .25(4?/ - 24 + .8y) = 16 - y' - fly. Rearranging, 
we find that 

4/ + 1.%' + fly = 16 

The initial conditions become y(0) = 0, and y'(0) = 6 + .25x(0) — .2y(0) = 6. 
Ordinarily as we'll see, it is more useful to convert higher order equations 
to lower. That always can be done. 

Problems 

1. Find a second order equation for xi, if 

x[ — X\ — X2 

x' 2 — Xl + 2x 2 

2. Find a second order equation for x 2 , for the same system as in problem 1. Is it the 
same as the equation for x\l 

3. Make a sketch of the phase plane for the system of problem 1. You should do this first 
by hand, and then let the computer do the work for you by running pplane on the system. 

4. Here you will find another system for the second order equation of problem 1. Set 
z\ = y, Z2 — y' — y, and find the equations for the derivatives of z\ and Z2- Did you get 
the same system as in problem 1? 

5. Make a sketch of the phase plane for the system 

x' = 1 

v = y - x 

and make a sketch of the slope field for the equation 

y = y-t 
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How do these two sketches compare, and why? 

6. Set up a system for a modified version of Jane's Candy Factory, in which there are now 
third and fourth tanks in series with the first two, with parameters as shown. 





tank 3 


tank 4 


sugar input 


lb/hr 


5 lb/hr 


water input 


15 gal/hr 





solution out 


25 gal/hr 


35 gal/hr 


tank contents 


220 gal 


165 gal 


weight of sugar 


z lb 


w lb 



You might notice that the volume of solution is no longer constant for all of the tanks. 

15 Systems of ODE's, part 2 

Today: Systems of differential equations, Newton, nonlinear springs, predator- 
prey, linear systems. Preview of algebra to come. 

There are several reasons for studying systems of differential equations. One 
is typified by Jane's Candy Factory of the previous lecture, where of course 
if we were accounting for 7 tanks, we would need 7 equations. Another case 
comes from Newton's F = ma law. Each mass m needs a second order 
equation, and there are many masses around. 

Example: Two Newton's Laws give 4 first order. Suppose we have two 
masses moving on the x axis pushing each other according to the system 

m 1 x" = f 1 (x 1 ,x 2 ) 
m 2 x 2 = f 2 {xi,x 2 ) 

In fact according to another of Newton's proclamations, f 2 = —fi, but we 
don't need this for our present purpose. Set vi = x[, v 2 = x' 2 . Then we have 
four variables representing the positions and velocities of the two masses. 
The system becomes 

Xi = Vl 
X 2 = V2 

miv'i = fi(xi,x 2 ) 
m 2 v' 2 = f 2 {xi,x 2 ) 



58 



Example: A nonlinear spring gives two first order. Consider the equation 
x" + 2x' + x 3 = 0. This equation looks somewhat like the second-order 
equations we studied for spring-mass systems, but it is not linear. You can 
think of the x 3 term as a nonlinear spring force. Specifically, in order to 
stretch this spring twice as far, you must pull eight times as hard. In many 
respects, this is closer to the behavior of a real spring. We now show how to 
convert this equation to a system of first order equations. Set y = x' . Then 
y' — x" = —2x' — x 3 = — 2y — x 3 . This gives us the first order system 

x = y 

V = -2y - x 

There must be some reason for doing this, right? First order systems do 
have some advantages. They are much easier to think about, for one thing: 
to imagine what happens with x" + 2x' + x 3 = involves thinking about 
the graph of x(t), and how the position x, slope x', and curvature x" all 
interact. This can be pretty complicated. On the other hand, thinking 
about an equivalent system just involves the idea that the solution curves 
must follow the vector field, as described previously. For emphasis, consider 
again a differential equation system 

x' = f(x,y,z) 

v = g(x,y,z) 

z' = h(x,y,z) 

The left sides give the velocity of a solution curve, while the right sides say 
that this velocity must agree with the vector field whose components are 
/, g, and h. In addition to this relatively simple description, most of the 
available software is also set up to deal with first order systems. 

Predator Prey Example Suppose we make a model of two populations, 
where x{t) and y(t) are the population sizes. We assume that these two 
species depend in each other in very different ways. As a first approximation, 
we assume that if there were none of the y species present, that the x species 
would grow exponentially, and that the exact opposite would happen to the 
y species if there were no x's around: we assume the y's would die out. Our 
first approximation is 

x = x first approximation : exponential growth 
y = — y first approximation : exponential decay 

This fits, if we are thinking of the y's as predators, and x's as prey, wolves 
and mice, for example. The next step is to add some terms for the interaction 
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of these species. We assume 

x = x — xy wolves are bad for mice 
y = — y + xy mice are good for wolves 

It should be clear that we have neglected to put in various coefficients; this 
does not affect the ideas involved. In a homework problem we ask you to 
run pplane on this system to see what happens. 

For now, we will just show that it is possible to begin with a system like this 
one, and convert it to a single second order equation, which is however quite 
a complicated equation sometimes. We calculate from the first equation, 
y = (x — x')/x = 1 — x~ 1 x'. This gives y 1 = x~ 2 (x') 2 — x~ 1 x' 1 ' . Then from 
the second equation, y' = (—1 + x)y = (—1 + x)(l — x~ l x'). The next step 
would be to equate these expressions and simplify. As you can see, this will 
give a very nasty second order equation. The first order system definitely 
makes more sense here. 

15.1 Linear Systems 

A linear system of differential equations means a system of the form 

x = 2x + 3y 
y = Ax + by 

The constants 2, 3, 4, and 5 might of course be replaced by any others. 
There could also be more than the two unknown functions x and y. You 
may notice that this seems to be a fairly simple system, since it does not 
include the predator-prey example, nor the nonlinear spring. We study 
these systems because 1) They are easy and fun to solve and make pretty 
pictures in the phase plane, 2) We learn some new algebra called linear 
algebra in the process of solving these, which is used for many things other 
than differential equations, and 3) They are the foundation for the harder 
equations. 

We'll attack the linear system above with exponential functions, since these 
have been useful so many times already. Try x = ae rt , y = be st . The system 
becomes 

are rt = 2ae rt + 3be st 
bse st = 4ae rt + 5be st 
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It is not quite clear what to do next, is it? But notice what happens if we 
take s = r. Then all the exponentials cancel 

ar = 2a + 36 
br = 4a + 56 

You may wish to compare this to what happened with second order linear 
differential equations. As in that case, we have gotten rid of the differential 
equations, and now have algebraic equations to solve. It should be easier 
this way. These are easier to solve if you rearrange them as 

(2-r)a + 36 = 
4a + (5-r)6 = 

What happens next can be confusing, if you don't see it coming. We are 
not going to solve this algebraic system right now. In fact, we are going to 
leave differential equations for a while to study the algebra connected with 
systems like these. It happens to be the same algebra which is used for 
many purposes not connected at all with differential equations, too. Once 
we know the algebra, we'll come back to the differential equation problem. 

Problems 

1. Find a first order system corresponding to x" — x + x = 0. You may use y = x' or 
something else that you choose, as the second function. 

2. Using y — x' , find a first order system corresponding to x" — x + x — 0. 

3. For the system you obtained in problem 2, substitute x = ae rt and y = be rt and find 
the algebraic system which results. You do not have to solve the algebraic system. 

4. Make a sketch of the phase plane for the system of problem 2. You should do this both 
by hand and by using pplane. 

5. Run pplane on first order systems for the equation x" + 2x + x p — for p = 1 and 
p = 3. What differences do you observe? 

6. Find a first order system corresponding to x" + x 3 — and try to sketch the phase 
plane by hand. This is too hard to do accurately, so don't worry about it. Next, multiply 
the equation by x' and integrate to find a conservation law, as in Lecture 10. Now it 
should be easy to sketch the phase plane, because the conservation law shows you what 
curves the solutions run along. 

7. Consider the predator-prey equations again. Our attempt to find a second order 
equation for this system was a mess, but it suggests something which actually works: you 
may notice the term ^ in part of our derivation. This suggests ln(x), doesn't it? since 
d(\n(x))/dt = Show that the equation x + y — ln(x) — ln(y) = c is a conservation law, 
in the sense that dc/dt — for all solutions to the predator-prey system. Sketch a graph 
of the function f(x, y) = x + y — ln(x) — ln(y). It may help to use octave for this sketch; 
type help or help plotxyz to find out how. 
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8. Here is way to "discover" the conservation law of problem 7. Imagine that the para- 
metric equations for (x(t), y(t)) are solved so that y is expressed as a function of x. Then 
by the chain rule you get ^ = ^j. This gives ^ = ~ y+xy , Show how to solve this by 

J J ° ax x ° ax x — xy J 

separation of variables. 

9. Run pplane on the predator-prey system and on the predator-prey system with mi- 
gration 

x = x — xy — .2 
y = -y + xy 

Which species is migrating here, and are they moving in or out? Write a verbal description 
of what the differences are for the two phase planes. Do you think there is a conservation 
law for the migratory case? 

10. What's rong with this? 

x = y 

l 3 

y = -x 

Then x" = -x :i , = -x 3 , ^ = dt 2 , % = dt + C, ln(x 2 ) = t, x = e^ 2 , y = 
(l/2)e t/2 + C. There are at least 3 mistakes. 



62 



16 Linear Algebra, part 1 



Today: Systems of linear algebraic equations. Matrices. 

A system of linear algebraic equations is something like 

2x - 3y + 5z = 
.7x — y + 3z = 
-y-2z = 7 

You might think, based on your past experience in math courses, that there 
is always exactly one answer to every problem. The situation here is a bit 
more subtle: This system can be thought of as three planes. If you plot 
these planes in an (x, y, z) coordinate system, you might find that all three 
happen to intersect in a single point. That point is the solution, in that 
case. 




Figure 16.1 Three planes on the left, and three on the right. When two 
planes in 3D intersect, they do so in a line unless they coincide. In one case, 
the three planes have only one point in common. In the other, something 
else happens. That is subtle in the picture. Can you see it? 

It may also happen that three planes meet in a line, rather than a point. 
Then there are infinitely many solutions. Again, it is possible that two of 
the planes are parallel. Then there is no solution. 



Practice: 1) There is another possibility. Can you think what it might be? 
2) Which of these possibilities holds for the system x + y = 0, 2x + 2y = 3, 

2 = 5? 
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16.1 Matrix 



A matrix is an array of numbers, such as the coefficients in our first system 



A 



2 
.7 




A matrix can also be formed to hold the variables 



u 



as well as the numbers on the right-hand side 

0" 

7 

Matrices are sometimes written using bold type, or with underlines, over- 
lines, overarrows, or rounded brackets, or square brackets around the name 
of the matrix. We won't use any of those. That means that when we use 
the symbol A or b, you have to keep straight from the context whether we 
are talking about a number or matrix or what. The size of a matrix is given 
as number of rows by number of columns. A is 3 by 3, while u and b are 3 
by 1. We introduce a new form of multiplication so that our system can be 
written 



or 



You might notice that to do this multiplication, it is helpful to move your 
left hand along the row of A, and your right hand down the column of u. 
The same concept is used to multiply matrices of any sizes, except that to 
multiply AB you must have the same number of columns in A as rows in 
B. Otherwise some poor number is left without anybody to multiply. 



"2 


-3 5 " 




X 




"0" 


.7 


-1 3 




y 










-1 -2 




z 




7 




Au -- 




b 







Example: A matrix multiplication: 



'2 3" 




X 




~2x + 3y 


4 5_ 




y. 




4x + 5y 
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1 2 

3 4 



2-3 5 
.7 -1 3 
-1 -2 

a b 
c d 



T 




"12" 







6.7 


2 




-4 



la + 2c 
3a + 4c 



16 + 2d' 

3b + 4d 



I is the special matrix 



one for each size. It is 



1 
or 1 or 
1_ 

the identity matrix and you should check that always AI = I A = A, Iu = u. 
It is also useful to abbreviate the system Au = b further by not writing the 
variables at all: the 3 by 4 matrix [A b] is called the augmented matrix of 
the system. In general, if you have a system of m equations in n unknowns, 
Au = b, then A is m by n, b is n by 1, u is m by 1, and the augmented 
matrix [A b] is m by (n + 1). 



16.2 Geometric aspects of matrices 



Before continuing with equation-solving, there are several remarks. Matrices 
like 

x 



u 



which consist of one column are sometimes called column vectors, and pic- 
tured as an arrow from the origin to the point with coordinates (x,y,z). 
In other words, u here has exactly the same meaning as xi + yj + zk, and 
you can think of it in exactly the same way. Of course we are now allowing 
our column vectors to contain more than three coordinates, so the three- 
dimensional vector notation doesn't carry over. Welcome to the fourth di- 
mension! There are some very nice geometric ideas connected to this form of 
multiplication. These seem at first very different and unrelated to equation- 
solving. We'll give one example now, just so you realise that there are several 
different reasons for studying matrices. Consider the simple matrix 



A 



'1 



A can be applied to any "point" (x, y) in the plane by using multiplication: 



A 



X 




X 


y. 




V 

.2. 
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This means that every point moves vertically closer to the x axis. In spite 
of how easy this multiplication is, the effect on the plane is dramatic, as you 
can see in the next figure. 




Figure 16.2 Matrices can change little smiley faces. What would happen if 
the face were not centered at the origin? 



Problems 
1. Let A = 



1 1 
1 



1 

1 



1 
1 
1 



D 



5 
5 



and let x = 



y - 



-l 

-2 

-4 



, B = V ' ( • 
Calculate AB, A 2 = A A, C 3 , DB, BD, Ax, 3Ax - 2Dz, Dx, and Cy. 



2. Suppose Au — b and Av = 6, and set s — u + v. Is it true that As = 6? How about if 
Av — instead, what then? 

3. An airline company plans to make x nights per month from its main hub airport to 
Chicago, and y flights to New York. New York flights require 6500 lbs of fuel and 200 
softdrinks, while Chicago flights require 5400 lbs of fuel and 180 softdrinks. The available 
monthly supply of materials at the hub consists of 260,000 lbs of fuel and 9,000 softdrinks, 
and all these are to be used if possible. Write a system of equations expressing all this. 
Discuss faults of this model. 



4. What's rong with this? 



"2" 




1 


2 




2 


4" 


4 




3 


4 




9 


12 
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17 Linear Algebra, part 2 



Today: Row Operations 

Now that our matrix notation is established, we turn to the question of 
solving the linear system Au = b. You may have solved such things in 
algebra courses in the past, by eliminating one variable at a time. This 
works fine for systems of 2 or maybe 3 equations. For our purposes though, 
it is important to learn a very systematic way of accomplishing this. We 
will describe a procedure which applies equally well to systems with any 
number of variables. In case you are wondering about the need for that, 
let me reassure you that linear algebra, of all the things you study, might 
be one thing that gets used in real life. Consider an airline company which 
needs to keep track of 45000 passengers, 200 planes, 300 crews, 8500 spare 
parts, 25 cities, 950000 pounds of jet fuel, 85600 soft drinks, etc, etc. Then 
there are the airframe engineers who compute the airflow at 10000 points 
around the plane, and the engine people who compute the temperature and 
gas flow at 3000 points inside the engines, and you get the picture. 

Before we describe row operations, you should see what the system of equa- 
tions looks like in a special case like 

"1 5" 
10 7 
1 9_ 

This is the augmented matrix for the system x = 5, y = 7, and z = 9. 
There is nothing to solve in this case. Consequently, the object of these row 
operations is to put the matrix into the form [/ b] if possible. 

We now describe "row operations" for solving the system Au = b. There 
are three row operations, which are applied to the augmented matrix [Ab], 
and have the effect of manipulating the various equations represented by the 
augmented matrix. 

1) A row may be multiplied by a non-zero number. 

2) Two rows may be interchanged. 

3) A row may be replaced by its sum with a multiple of another row. 

Row operations do not change the solutions of the system. This is clear for 
the first two, and requires a bit of thought for the third one. They only 
change the appearance of the equations. 
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A system may have many solutions. 



"2 


3" 




X 




r 


4 


6 




y. 




2 



Example: The system 



represents the two lines 2x + 3y = 1 and 4x + 6y — 2. You might notice that 
these are two equations for the same line. Therefore any point on this line 
should be a solution to both equations. Let's do the row operations and see 
what happens. We have the augmented matrix 



2 3 1 
4 6 2 



(row 2 — 2 row 1) 



2 3 1 




Here we have replaced row 2 by (row 2)-2(row 1). Look at the bottom line. 
It says that Ox + Oy = 0. That is true but not helpful. The top line says 
that 2x + 3y — 1 . Since this is the only requirement for a solution, there are 
indeed infinitely many solutions 



2 2 y 
y = anything 



A system may have no solutions. 



Example: The system 



has augmented matrix 

2 3 
4 6 2 



'2 3" 




X 




"0" 


4 6 




y. 




2 



(row 2 — 2rowl) 



2 3 
2 



Look at the bottom line. It says that Ox + Oy — 2. Is that possible? (No.) 
So this system has no solution. It describes two parallel lines. 

The next example involves a bit more arithmetic, and is our original system. 
It shows that a system may have a unique solution. 



Example: 





"2 -3 5 


0' 




"1 -1.5 2.5 0" 




[Ab} = 


.7 -1 3 





— > (.5 row 1) — > 


.7-1 3 


— > (row 2 




_0 -1 -2 


7 




-1-2 7 






"1 - 


-1.5 2.5 0" 















05 1.25 


— > (interchange row 2 and row 3) and ( 


—row 2) — > 







-1 -2 7 
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1 -1.5 2.5 
1 2-7 
.05 1.25 



1 5.5 
1 2 
1 



1 5.5 
12 
1.15 

-10.5 ^ 
-7 
.3043 



(row 1 + 1.5 row 2) and (row 3 — .05 row 2) 
(row3)/1.15 -> 



-10.5 
-7 
.35 



(row 2 — 2 row 3) and (row 1 — 5.5 row 3) 



1 -12.1737 
1 -7.6086 
1 .3043 

The solution to our original system is therefore 

-12.1737 
-7.6086 
.3043 

Notice there is only one solution, which says that these three planes intersect 
in just one point. You may notice that this is the system graphed in Figure 
17.1, on the left. 



A matrix is "row-reduced" if it looks like one of 



1 
1 



"1 





0" 




"1 


4 





5" 





1 













1 


7 








1 

















1 3.55 
1 2.17 



, etc. 



i.e., the leading non-zero in each row is 1, the column of that 1 contains 
no other non-zeroes, and these leading l's move to the right as you read 
down successive rows of the matrix. Any row of zeros is at the bottom. A 
row-reduced augmented matrix is easy to solve: For example the 3 by 4 
matrix represents the system 

x + 4y = 5 
z = 7 
= 



which has solutions 



7 

anything 
5-4y 
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The downright pattern of l's allows you to solve the reduced system from 
the bottom up. 

We point out again the term "linear" . Multiplication of vectors by a matrix 
is linear in the sense that 



A{c\ui + c 2 n 2 ) = c\Au\ + c 2 Au 2 

where the c's are any constants. Consequently if you have two solutions 
u\ and U2 to Au = 0, it follows that some other solutions are the so- 
called "linear combinations" c\U\ + c 2 u 2 . So we say that Au = is a linear 
equation, and this is linearity in the strictest sense of the word. We also say 
that Au = b is a linear equation, but here you have to be careful to see that 
a linear combination of solutions is usually not a solution. 

Problems 

1. Let A = q ^ and B = J j . What single row operation has been done to A to 
get B? 

2. What row operation can be done to A of problem 1, to produce the identity matrix 7? 



3. The matrix 



needs two more row operations done, to put it into the 



10 3 7 
2 1 

reduced form. What are they, and what is the reduced form? Note that the reduced form 
is unique, but there are two different sets of row operations which are capable of producing 
it. 

4. Write out the system of algebraic equations corresponding to the augmented matrix 



Is there something odd about the third equation? How many solutions does this system 
have? 

5. Write the augmented matrix for the following system, row reduce it, and find all the 
solutions. 

x — y — 

x + y = 1 

—x — 3y — —2 



6. The following is a reduced augmented matrix for a system of 4 equations in the five 
unknowns x\ . . . x 5 . Find the solutions to the system. 



1 





3 














1 











4 











1 





5 














1 
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7. Same instructions as problem 5 for the system 



x — y + z — 2w = 
x + y + z — w — 
-x-y = 




1 
2 



8. Whats rong with this? 



x + 2y = 
3y = 



gives y = 0, then x = 0. Therefore there are no solutions. 

18 Linear Algebra, part 3 

Today: Matrix inverse, solutions to linear sytems in matlab, and the Eigen- 
value concept 

18.1 Matrix Inverse 

We learned last time about solving a system Ax = b by row operations - you 
row reduce the augmented matrix [A b] and read the solutions, if there are 
any. There is another approach which is harder to compute but helpful to 
know about sometimes. It goes like this. We try to solve Ax = b by analogy 
with the arithmetic problem 2x = 7. In arithmetic the answer is x = 
2~ 1 7. Here there is sometimes a matrix which deserves to be named A^ 1 , so 
that the solution might be expressed as x = A~ 1 b. A" 1 is pronounced "A 
inverse". Wouldn't that be nice? Yes, it would be nice, but we'll see some 
reasons why it can't work in every case. It can't work in every case because, 
among other things, it would mean that there is a unique answer x, and we 
have already seen that this is not at all true! Recheck the examples in the 
previous sections if you need a reminder of this point. Some systems have 
no solutions, some have one, and some have infinitely many. It turns out 
that for some matrices A, the inverse A -1 does exist. All of these matrices 
are square, and satisfy some additional restrictions which we will list below. 
Before getting to that list though, return to the phrase used above, "a matrix 
which deserves to be named A" 1 " . What does that mean? We'll say that 
a matrix B is an inverse of matrix A if AB = BA = /. This is the key 
property. In fact, there can only be one matrix with this property: if B and 
C both have this property, then B = BI = BAC = IC = C. Since there is 
only one, we can name it A -1 . 
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Example: For A 
and BA = I. So 

and 



1 



and B 



1 
2 



1 

2 



1 
2 



1 

i 



we have AB 



1 
1 



Example: For A = 
I and BA = I. So 

and 



and B = 



2 5 
1 3 

' 2 
-1 



3 

-1 



, we have AB = 



2 5 
1 3 



1 
1 



It turns out that there is a formula for the inverse of a 2 by 2 matrix, but 
that there is no really good formula for the larger ones. It is: 



a b 
c d 



d 

— c 



ad — be 

You can derive this formula yourself: try to solve 



IF ad-bc^ 0. 



a b 




x y 


c d 




z w 



= I for 



x, y, z, and w in terms of a, b, c, and d, and you will get it eventually. 
The number ad — be is called the "determinant" of the matrix, because it 
determines whether or not the matrix has an inverse. 



Practice: You should check that this formula for the inverse agrees with 
the earlier examples. 

We also are going to record the method used to find inverses of larger matri- 
ces. They may be found, when they exist, by row operations. You write the 
large augmented matrix [A I] which is n by 2n if A is n by n. Then do row 
operations, attempting to put it into the form [IB]. If this is possible, then 
B = A^ 1 . If this is not possible, then A' 1 does not exist. The rationale for 
this method is left for a later course in linear algebra, but you can consider 
for yourself what equations are really being solved, when you use this large 
augmented matrix. 
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Example: For A = 



-10 
2 1 
4 



[A I] 



-10 10 
2 10 10 
4 1 



we do operations: 



— (rowl) and (l/4)(row3) and (row 2 + 2rowl) 



10 0-10 
10 2 10 
1 .25 





'-1 





0" 




"-1 





" 


You can check that 


2 


1 







2 


1 













4 










.25 



= / 



Have you heard about people who "only know enough to be dangerous" ? It 
means they have some special knowledge of something, but without enough 
depth to use it properly octave and matlab provide a tool to put you into 
this category Isn't that exciting? So we hope that you will not get hurt 
if you know the following: the program can sometimes solve Au = b if you 
simply type u = A\b. That is a backward division symbol, intended to 
suggest that one is more-or-less dividing by A on the left. It mimics the 
notation u = A~ l b. In the problems, we suggest trying this to see what it 
produces in several cases. 



18.2 Eigenvalues 

Usually when you multiply a vector by a matrix, the product has a new 
length and direction, as we saw with the little smiley faces: If u is a vector 
from the origin to the left eye on page 66 then Au points in a different 
direction from u. This is typical. On the other hand, let v be the vector 
from the origin to the top of Smiley 's head. In this case, we see that Av = \v. 
This is not typical, it is special, and turns out to be quite important. Try 
to locate other vectors with this special property, that 

Av = Xv 

A vector v ^ which has this property, for some number A ("lambda"), is 
called an eigenvector of A, and A is called an eigenvalue. Eigen is a German 
word which means here that these things are characteristic property of the 
matrix. You may notice that we excluded the zero vector. That is because 
it would work for every matrix, and therefore not be of any interest. It is ok 



73 



for the eigenvalue to be 0, though. That just means that Av = for some 
nonzero v. Since this is not true for all matrices, it is of interest. 

Note that this eigenvector equation Av = Xv is superficially like the system 
equation Ax = b, but there is a big difference. The difference is that b is 
known, whereas the right-hand side of the eigenvector equation contains 
the unknown v and the unknown A. So the eigenvector equation is probably 
harder. Let's write out the eigenvector equation for the matrix 



A 



2 3 
4 5 



If we write a and b for the components of v, then Av = Xv becomes 



"2 


3" 




a 


= X 


a 












4 


5_ 




b_ 


b 



or 

2a + 36 = Xa 
4a + 5b = Xb 

Now go back to reread the discussion of the linear system of differential 
equations near the end of Lecture 15. Do you recognize the equations? 



Example: Reconsider the matrix A — 



1 

as on page 66. Any vector along the x\ axis is unchanged, and is therefore an 



eigenvector with eigenvalue 1: A 
vectors? 



which squeezes smiley faces 
langed, and is therefore an 
What happens to vertical 



Problems 
1. Solve 



three ways: by row operations, matrix inverse, and the Matlab \ command. 
2. Same as problem 1, for 



2 


3" 




6 


4 


5_ 


u = 


7 



"2 


3 





0" 




"0" 


4 


5 








u — 











1 


1 





4 


6 


1 







1 
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3. Find the inverses of the matrices 



1 a 
2 



and 






1 














1 














1 


a 












Is there any restriction on a in either case? 





"1 


2 


0" 




"1" 




"3" 




"0" 




"0" 




'l" 


4. Given A = 


2 


4 





, u = 


2 


, v = 


6 


, w = 





, x = 





, y = 


l 










5 














1 




7 




l 



and 



-1 







Which of u, v, w, x, y, and z are eigenvectors of A, and what is the eigenvalue 
in each case? 

5. Find all the eigenvalues and eigenvectors for all the identity matrices I. 

6. Run matlab on the three examples of Lecture 17: to solve Ax = b, set x = A\b. 
Warning: We know ahead of time that this must fail in at least two of the three cases. 
Discuss the results. 

7. Write down a matrix which stretches smiley faces to 4 times their original height. Find 
all the eigenvectors and eigenvalues of this matrix. 

■i i" 



8. Whats rong with this? 



b, x = 



19 Linear Algebra, part 4 

Today: Eigenvectors and Eigenvalues 

We have seen the eigenvector concept: a non-zero vector x is an eigenvector 
for a matrix A with eigenvalue A if Ax = Ax, and we have seen some exam- 
ples. In this lecture we will outline some methods for finding eigenvectors. 
For our purposes in this course, we need to know the 2 by 2 case very well, 
the 3 by 3 somewhat, and the larger sizes just a little. We begin with a 2 

by 2 matrix A = a ^ . Consider the following statements. 

1. A does not have an inverse. 

2. det A = 0. 

3. One row of A is a multiple of the other. 

4. is an eigenvalue of A, i.e. there is a non-zero x with Ax = 0. 
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5. Solutions to Ax = b, if they exist for a particular b, are not unique. 

These statements are all equivalent, in that if one of them is true of a partic- 
ular matrix A, then they all are true of A. We discussed the equivalence of 1 
and 2 on page 72. We won't prove all the equivalences, but the implication 
(4 implies 1) is the most interesting for our purpose, so let's think about 
that one: assume that there is a non-zero vector x such that Ax = 0. Now 
the claim is that A cannot have an inverse under this circumstance. The 
proof is very easy. Suppose A^ 1 exists. Multiply the equation Ax = by 
A^ 1 to get A -1 Ax = A^ 1 or x = 0. This contradicts the fact that x is 
non-zero. Therefore A~ l does not exist. 

Now we pull ourselves up by our own bootstraps, using these facts about 
eigenvalue to help us find all eigenvalues. So suppose now that B is any 2 
by 2 matrix and we want to solve 

Bx = Ax 



Note first that this may be rewritten as 

(B - XI)x = 

where I is the identity matrix. Now we just apply the previous discussion, 
taking A to be (B — XI). The conclusion is that 

det(JB - XI) = 

This is our equation for A. You can see, or will see, that it is a quadratic 
equation. Once it is solved, we solve the previous one for x. 



Example: 



det(B - AJ) = det( 



B = 



1 2 

-3 



A 



det( 



1 - A 




2 

-3- A 



(1 



A) (—3 — A) — = 0. This is already factored so we find A = 1 or A = —3. 
These must be the eigenvalues of B. The next step is to find the eigenvectors. 
Take the case A = 1. We must solve (B — ll)x — 0. Explicitly this says 



"0 


2 " 




Xl 










-4 




X2_ 








These require X2 = 0, and note that is the only requirement. So you can take 

Xl 







where xi is anything other than 0. Usually we take the simplest 
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possible thing, which is 



. We have found, tentatively, 



B 



1— l 




1 











and you should check this. Then do the case A = —3: We must solve 
(B + 3I)y = or 



We 



4 


2 




yi 




"0" 










2/2 








These require j/2 = — 2j/i. So you can make a simple choice y 
have found 

B 



Again, you must check this. Does it work when you multiply it out? Ok. 
Note that we took simple cases, but that there are also other choices of 



eigenvectors which are multiples of these, such as 
Check these too if you're not sure. 



5 

-10 



and 



-3.172 
6.344 



Note that you can make a picture of the plane showing all the eigenvectors 
of B, and you will find that they lie on the two lines y = and y = —2x, 
excluding 0. These are called the eigenspaces for the eigenvalues 1 and —3, 
respectively, except that the vector is included in the eigenspaces even 
though it is not an eigenvector. Everything on the x axis remains unmoved 
by B, while everything on the line y = —2x gets stretched by 3 and flipped to 
the other side of the origin. Note that having all this eigenvector information 
gives us a very complete understanding of what this matrix does to the plane: 



Practice: What would happen to a smiley face under the action of B1 



Example: 



C = 



1 

-9 



We get det(C - AI) = A 2 + 9. This is zero only if A = ±3i. For the case 
A = 3i we solve (C — 3il)x — which reads 



= 



~-3i 


1 




Xl 


-9 


-3i 




X2_ 



So X2 — 3ixi . We take x 



Then check 



C 



"1" 




"1" 




= 3i 


3i 


3i 
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For the case A = — 3i, the only change is to replace i by —i everywhere, and 
check that 



C 



1 " 




1 " 




= -3i 




-3i 


-3z 



Note that for C, one has run into complex numbers, and so the nice ge- 
ometric interpretation which was available in the plane for B is no longer 
available. Two complex numbers are about the same as four real numbers, 
so you might have to imagine something in four dimensions in order to 
visualize the action of C in the same sense as we were able to visualize B. 



Practice: It does make sense, nevertheless, that not all 2 by 2 matrices 
can be thought of as stretching vectors in the plane: figure out what what 
C does to smiley faces, if you want to understand this point better. 



Larger Matrices 

The procedure for 3 by 3 and larger matrices is similar to what we just did 
for 2 by 2's. The main thing which needs to be filled in is the definition of 
determinant for these larger matrices. We do not want to get involved with 
big determinants in this class. To find eigenvalues for a 3 by 3 matrix you 
can look up determinants in your calculus book, or remember them from 
high school algebra. Also see Lecture 20. There is also a possibility of using 
some routines which are built into octave. These work as follows. Suppose 
we have a matrix 

"1 2 3" 

A= 4 5 6 
_7 8 9_ 

To approximate the eigenvalues of A you can type the following commands: 

A=[l 2 3;4 5 6;7 8 9] ; 

eig(A) 

To get eigenvector approximations also you may type 
[V D]=eig(A) 

This creates a matrix V whose columns are approximate eigenvectors of A, 
and a matrix D which is a "diagonal" matrix having approximate eigenvalues 
along the main diagonal. For example, D in this case works out to be 



D 



16.1168 







-1.1168 
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Problems 



1. Sketch a smiley face centered at the origin and its image under the matrix B of the 
text. 

2. Repeat problem 1 for the matrix C of the text. 

3. The last example 



A = 



1 2 3 
4 5 6 
7 8 9 



seems to have as an eigenvalue. Find a corresponding eigenvector. 
4. Find eigenvalues and eigenvectors for the matrix X = 



2 -3 

3 2 



5. Find eigenvalues and eigenvectors for the matrix W = 



1 

-1 -1 



6. Which of the vectors 



1 




1 




l 




1 


2 




2 




-2 




2 


2 


1 


1 




-2 




2 


1 




2 




1 




-1 



are eigenvectors for the matrix 









1 














1 


1 














1 









and what are the eigenvalues? 

7. Which of the vectors in problem 6 are eigenvectors for /, the 4 by 4 identity matrix 
and what are the eigenvalues? 

8. Which of the vectors in problem 6 are eigenvectors for 77 and what are the eigenvalues? 

9. Which of the vectors in problem 6 are eigenvectors for the matrix 



1 








0" 








1 








1 

















1 



and what are the eigenvalues? 

10. If you don't know how to do the determinant of a 3 by 3 matrix, look it up somewhere, 
or see page 80 and find the eigenvalues and eigenvectors for 



A = 



3 

5 



3 
5 7 



11. What's rong with this? A = 
A) = 21 - 7A - 3A + A 2 = A 2 - 10A + 21, so use the quadratic formula 



det(A-XI) = det( 



3- A 
5 





7- A 



(3-A)(7- 
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20 More on the Determinant 



Today: The determinant determines whether the matrix is invertible. 



Unlike most of the lectures we aren't going to explain everything here, but 
will indicate most critical ideas. 

2x2 

You know that for 2 by 2 matrices, the determinant is 



det 



a b 
c d 



ad — be 



and it determines whether the matrix is invertible. It is used in the formula 
for the inverse 



a b 
c d 



-i 



ad — be 



d -b 
— c a 



It also tells geometric effects, for example if det(A) = —3 then A expands 
areas by a factor of 3 and reverses the orientation. The determinant was 
discovered when someone (like you, you could try this) solved a linear system 
symbolically and realized that (ad— be) always shows up in the denominator. 

3x3 



If you solve the system 



a b c 
u v w 



r s t 

symbolically then you find the expression 



avt + bwr + cus — aws — but — cvr 

in the denominators of the x^. That expression is then defined to be the 
determinant of the 3 by 3 matrix. It is so complicated that people have 
found various ways to express it, such as: 



det 



a b c 
u v w 
r s t 



a det 



v w 
s t 



b det 



u w 
r t 



+ c det 



u v 
r s 
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Example: det 



= 3. This matrix stretches or compresses the 



3 
.5 
2_ 

three axes by various factors, and it multiplies volumes by 3. The matrix is 
invertible. The determinant of the inverse is |. 

Any Size 

You can find the determinant of any square matrix iteratively as in this 
example, note the alternating signs: 



det 



a bed 

12 

-2 4 6 7 

-3 5 8 



a det 





2 





4 


6 


7 





5 


8 





" 


2 


0" 




" 


1 


0" 




" 


1 


2" 


b det 


-2 


6 


7 


+ c det 


-2 


4 


7 


— d det 


-2 


4 


6 




-3 


5 


8 




-3 





8 




-3 





5 



Then do the 3 by 3's, etc, like: 



det 



1 2 
4 6 7 
5 8 



= 1 det 



6 7 

5 8 



- 2 det 



4 7 
8 



det 



= 1(48 - 35) - 2(32 - 0) + 0(20 - 0) = -51 



20.1 Products 

Another generally useful fact: &et(AB) = det(A) det(B). For example, if A 
is invertible then 

1 = det(J) = detiAA- 1 ) = det(A) det(A- 1 ) 

so det(A- 1 ) = gA^. 

To prove the product rule is for another course. But we can use it to prove 
something good: 

Cramer's Rule: M. Cramer cleverly noticed for example that 



If 



a b c 
u v w 
r s t 



Xl 




~P 


X2 




D 


X3_ 




Q 



then 



a b c 
u v w 
r s t 



xi 0' 
x 2 1 
x 3 1 



P b c 
D v w 
Q s t 
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and taking determinants you get a formula for the answer: 

P b c 
D v w 
Q s t_ 

Similary for X2 and xs, and for any size system. Cramer's rule is an inter- 
esting theoretical tool these days, but not computationally efficient because 
it takes too long to work out determinants. 

Problems 

1. Is this matrix invertible? 

"0 8" 
10 
10 
1 0_ 

2. If a 2 by 2 matrix A has columns v and w we will write A = [v w] and think of the 
determinant as a function of the two columns f(v,w) = det(A). Convince yourself that 

f(au + bv, w) = af(u, w) + bf(v, w) 

whenever u, v, and w are vectors and a and b are numbers. 

We say that det is a linear function of the first column of the matrix. 

3. Of course det is a linear function of the second column too. Figure out why that makes 
this kind of identity hold: 

f(au + bv, cw + dz) = acf(u, w) + adf(u, z) + bcf(v, w) + bdf(v, z) 



det(A)xi = det 
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21 Linear Differential Equation Systems 

Today: Just when you thought it was safe! back to linear differential equa- 
tions. The solution by eigenvectors. 

At this point you have seen those essential parts of linear algebra which are 
needed for solving linear systems of differential equations of the form 



Here, A may be an n by n matrix and x an n by 1 column vector function of 
t. As in any differential equation system x' = f(x), the right side represents 
a vector field which exists in space, and the mission of any curve x{t) which 
intends to solve the system is to be sure that it's tangent velocity vector 
x'{t) agrees at all times with the given vector field. In this lecture we deal 
only with vector fields of the form f(x) = Ax, and these may be solved 
rather completely using our new eigenvector knowledge. 



Figure 21.1 A vector field in 3 dimensions. We think of it as unchanging 
unless we have a need for explicit time dependence in the field. But there 
are also two curves shown, which you might think of dynamically. They are 
doing their very best to follow the arrows, not only in direction but also to 
have the correct speed at every point they go through. If they succeed in 
this, they are solutions to an ODE. The field shown is more complicated 
than a linear Ax field. 

For example, we know that sometimes you find an eigenvector v for which 
Av is a positive multiple of v. In such a case, the vector field points away 



x' = Ax 
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from the origin and so any solution passing through point v must be headed 
away from the origin. Similarly negative eigenvalues give rise to solutions 
headed for the origin. The effects of complex eigenvalues are a little harder 
to guess. So let's proceed with the analysis in the following way. Suppose 
by analogy with many things which have gone before, that we are looking 
for a solution of the form x(t) = ce rt , where r is a number as usual, but now 
c must be a vector for the assumption to make any sense. Then substitute 
into the differential equation to get 

x' = rce rt = Ax = Ace rt 

These will agree if 

Ac = rc 

Now stop and take a deep breath, and think about what you learned before 
about linear algebra. This equation is the same except for notation as 

Av = Xv 

isn't it? We have reached the following fact: 

If Av = Xv, then one solution to x' = Ax is ve xt . 



Example: From the example on page 76, 
solution to 

x — x + 2 
V = Sy 



'1 2 " 




T 


= 1 


"1" 


-3 













So one 



e 1 * i.e., x = e*, y = 0. Check these in the differential equation to be 



Note that we certainly don't stop with one solution. For each eigenvector 
you get a solution. 



Example: Also for the same example we found that 



is an eigenvector 



with eigenvalue —3. Consequently we have another solution 
x = e~ 3t , y = -2e~ 3t . Check these too. 



1 

-2 



In general we hope for a system x' = Ax where A is n by n, that there 
will be n such solutions. Usually there are, but not always. More on this 
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in a moment. The most important thing now is to point out how you can 
combine these solutions which we have found. Our equation x' = Ax is 
linear in the strictest sense of the word: linear combinations of solutions are 
solutions, i.e., if x[ = Ax\ and x' 2 = Ax2, then for any constants c\ and C2 it 
is true that (c\X\ + C2X2)' = A{c\X\ + C2X2). This is because differentiation 
and multiplication by A are both linear operations. You should write out 
the one-line proof of this yourself if you are not sure about it. 



Example: x(t) 





T 




" 1 " 


Cl 





e* + c 2 


-2 



is a solution to x' 



1 2 

-3 



for any choice of constants ci and C2 . In order to meet a given initial condi- 

[2.5" 

tion, we have to solve for the Cj. For example if we want to have x(0) = 



2.5 




1 


+ C 2 


1 


-2 


= Cl 





-2 



then it is necessary to solve 
ci = 1.5 and we get the solution x(t) = 1 



The answer is 01 = 1, 



e 1 + 1.5 



Back to the question of finding n solutions for the case of n by n matrices 
A. Certainly if we find n eigenvectors Vj we can form solutions like 



x(t) = ci«ie Al * + 



The question is whether these are all of the solutions to the differential equa- 
tion. Equivalently, is it possible to meet all the possible initial conditions 
using 

ar(0) = C1V1 -\ h CnV n ? 

Note that this is a system of linear algebraic equations in the unknown Cj, 
like we saw in the previous sample. What is needed here is that the Vj 
should be linearly independent. This concept is slightly beyond the level of 
these notes, but basically it means that the vectors really point in n different 
directions, or that none of them may be written as a linear combination of 
the other ones. You will study this concept in a linear algebra course; we're 
only doing an introduction to it here. 



Example: Let A 



3 

5 



This was problem 10 on page 79. We have 





T 




"0" 


Avi — 3vi and Av2 = where vi = 





, V 2 = 












1 



, and there are no other 
eigenvectors except for multiples of these. So some solutions to x' = Ax are 
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x — ciVie M + C2V2. This is one of the exceptional cases where you do not 
get n(= 3) eigenvectors, nor the general solution to the differential equation. 
See if you can find another form of solution to these differential equations; 
for this system it is actually easier to write the equations out and deal with 
them explicitly, rather than to think about eigenvectors. 



Problems 
1. Solve 

x = 2x 

V = % 

x(0) = 5 

2/(0) = 



2. Solve the system of problem 1 with the initial conditions 

x(0) = 
1/(0) = 5 



3. Sketch the phase plane for the system of problem 1. 

4. Solve the system of problem 1 with initial conditions 



x(0) = 
2/(0) = 



Sketch your solution onto the phase plane, or run pplane on it. Show by eliminating t 
between your formulas for x(t) and y(t) that x = j/S. This should help your sketch. 



5. Solve 



2x 



-32/ 



Show that solutions lie on curves x — cy 3 , or lie along the positive or negative axes, or 
are zero. Sketch the phase plane. 



6. Sketch the phase plane for the system 



which was solved in the text. Does it look more like the phase plane for problem 3 or 
problem 5, qualitatively? 



7. Solve 



and sketch the phase plane. 
8. Solve 



x — 
i 

y = 



x + 2y 
Ax — y 



86 



y = x + 5y 

and sketch the phase plane. 

9. You are given that a matrix A has eigenvectors u 



1 















, v = 


3 


, w = 


2 







-1 




-2 



and that Au = 5u, Av — 4v, and Aw — 3w. Solve x' = Ax, if x(0) 

x(0) =u + 2v. 

10. What's rong with this? 



and A 



T 




'-7' 


2 




-14 



1(0) = 
is given. Soi = 3e~ 7t + C. 



Ax 
3 
6 



Note that 



22 Systems with Complex Eigenvalues 

Today: A complex example. 

Today's lecture consists entirely of one example having complex eigenvalues. 
It is not really that hard, but we want to be careful and complete, the first 
time we work through one of these. Our system is 



where 



= Ax 

1 2' 

-2 1 



We know what has to be done to solve this; you find eigenvalues then eigen- 
vectors of A, then the solution x(t) is in the form of linear combinations of 
the ve xt things. This is what we have done before, and here we will find a 
bit of complex arithmetic to go along with it. Here we go. 

First plot the phase plane to see what to expect. You can plot a little bit 
of it by hand, or you can run pplane on our system: 

V = V + 2z 
z = -2y + z 



87 



2- T 

t 

T 



t r 

T T 

T T 

r t 

1 T 



? 7 ^ 



' *- *r- ^ 

^ ^ ^ 
~ ^ ^ ^ ^ ^ 
• * . * . 



* -* ^* -* ^ iSt ^ Si 

«»~. A, ^IVi \ 

» -» "» V "a ^ \i \, \i \ 

•^%.^.-*\\\\\\ 

• ^ V \ ^, \, \, \ V \, \ - 

y ^ i j j j i i j i - 
^///^ J Jiii 

<s / ^ ,/ j. ,i ^ ,i j, 
s ^ 4 ./ <l J, .1 J„ -I - 



Figure 22.1 The phase plane for our example, (figure made by pplane) 



Now let the analysis begin. First we find the eigenvalues of A. We have 



det(^ - XI) = det( 



1 - A 2 
-2 1 - A 



) = A 2 - 2A + 5 



This is zero by the quadratic formula if A = 2± ^\ = ~ 
So we have our eigenvalues. 



-16 



= 1 ± 2i 



Now let's look for an eigenvector for l + 2i. We must solve (^4 — (1 + 2i)I)v = 
0. Writing a and b for the components of v, this says 

-2ia + 26 = 
-2a -2ib = 

The first equation says b = ia. The simplest choice we can make is to take 
a = 1, b = i. Now something strange happens. If you check these numbers 
in the second equation, you will find that they work there too! In fact, you 
can just forget about the second equation, because it is a multiple of the 
first one. I know it might not look like a multiple of the first one, but that is 
because it is an imaginary multiple of the first one. Anyway we have found 
that 



1 

-2 



(1 + 2*) 



Read that through again to check it, because there are lots of places you 
can make mistakes in complex arithmetic. 
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Next let us state the solution we have found so far, and expand it so we can 
see what it looks like. 



e (1+2i)t = e*(cos(2t) + isin(2t)) 



cos(2t) ' 
- sin(2t) 



le 



sin(2t) 
cos(2i) 



That's a mess, isn't it? It is neater to leave it in the original complex 
exponential form, but we expanded it out so you can see exactly what the 
real and imaginary parts are. Now we need a second solution, so we can go 
through all this again with the other eigenvalue l — 2i. The only thing which 
changes, if you read through what we just did, is that you must replace i 
by —i wherever it occurs. This happens to work because A is a real matrix. 
The general solution to the system can then be expressed as 



x(t) = Cl 



e( 1+2i >* + C2 



1 



3 (l-2i)t 



There is an alternate to this, sometimes preferred because it emphasises 
real solutions. This goes back to our discussion in Lecture 12 about real and 
imaginary parts of solutions to real equations. The same statement holds 
here, that if x = u + vi is a solution to x' = Ax where u, v, and A are real, 
then u and v are also solutions. You can prove this yourself pretty easily. 
Taking the real and imaginary parts we found for the first solution above 
allows us to express solutions in the alternate form 



x(t) = a\e l 



cos(2t) 
- sin(2t) 



+ a 2 e l 



sin(2i) 
cos(2t) 



The relations between the constants are a± = c± + c 2 and a 2 = i(c\ — C2). 

It is important to think about the phase plane and see whether it matches 
the solutions we found. You can see that we got the increasing exponential 
e*, and this has to do with the fact that solutions are moving away from the 
origin in Figure 22.1. Also we got sines and cosines, which oscillate positive 
and negative. That has to do with the rotational motion in the picture. 

Problems 



1. Solve x' = Ax if A 



-2 
1 



Find the solution if x(0) = 



2. Solve x' = 
in this case? 



1 

-9 



x. What is the second order equation for the first component xi 
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3. Write out the system corresponding to x" = —x using x' = y and solve it, find a 
conservation law for the equation and plot its level curves, and plot the phase plane for 
the system. How are these plots related? 

4. Solve 



x 

t 

y 



= —3a; 

= 2 

= -9y 



and plot some solutions. 



with this? x — 


T 


e (l+2i)t = 


"1" 




i 




i 



e*(cos(2t) + isin(2t)), so re(ai) 



e* cos(2t), im(x) = e* sin(2t), and x = e'(ci cos(2t) + C2 sin(2t)). 
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23 Classification Theorem for Linear Plane Sys- 
tems 



Today: Classification of linear planar systems. Spirals, centers, saddles, 
nodes, and other things more rare. 

Today we have a nice treat for you. Rather than doing a lot of computa- 
tion, which we have been doing for several lectures, we are going to state 
a classification theorem. If this sounds like when the doctor says "It won't 
hurt at all", that's because we first have to explain the significance of what 
we are about to do. A classification theorem is one step higher in the world 
than computation, because it is something which lists all the possible things 
which could ever happen, in a given circumstance. Having such a theo- 
rem is considered by mathematicians to be a good thing. There are not 
very many classification theorems, and mathematicians wish they could find 
more. Here is what this one says. 

Theorem on Classification of Linear Systems in the 
Plane Suppose you have a system x' = Ax where A is a 2 by 
2 real matrix. Draw the phase plane. Then the picture you just 
drew has to be essentially one of the seven pictures shown below, 
and there are no other possibilities. Essentially means that you 
might have to rotate or stretch your picture a little to make it fit 
or reverse the arrows, but there will be no more difference than 
that. 

There now, that was not so bad, was it? Here are the pictures. The most 
important of these should look familiar from homework. The classification 
proceeds by considering what the various possibilities are for the eigenvalues 
of A. 
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Figure 23.1 This is called a "node". It is what happens when the eigen- 
values are both real numbers, and have the same sign. The case shown has 
both eigenvalues positive. For both negative, the arrows would be reversed. 
Important: solutions do not go through the origin. 




Figure 23.2 This is called a "saddle" , because it is the pattern of rain running 
off a saddle, as seen from above. It occurs when the eigenvalues are real, but 
of opposite sign. 
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Figure 23.3 This is called a "center". It occurs when the eigenvalues are 
±bi, for some non-zero real number b. 
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Figure 23.4 This is a "spiral", for obvious reasons. The eigenvalues are a±bi 
with a 7^ 0. The case a > is shown, but if a < the arrows are reversed. 
Also this can appear clockwise or counterclockwise. 



93 
















1.5 


























< 










1 
























0.5 


- it? 

























_ ^ 


1^ 














\S 
















0.5 
























-1 













Z y? Z i? Z 

5 L*r' ^ <Z ^ <S Z ^ 

<S S I? ^Z 

~ 2 r< ^ iz ^ ^ ^ 

-2 -1.5 



Figure 23.5 A rare case in which one eigenvalue is 0. 
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Figure 23.6 A rarer case in which both eigenvalues are 0. 



94 



Figure 23.7 The rarest case, in which the entire matrix is 0. 
Problems 

1. If you did problems 1,2, and 3 of the previous lecture, determine which of our phase 
planes in this lecture applies to each of them. 

2. Make up your own system x' — ax + by, y' = cx + dy, by choosing the coefficients 
at random. Draw the phase plane for your system by hand or pplane and classify which 
type it is by trying to match it with one of the figures in this lecture. The theorem says 
that it should match. 

3. Is there any other classification theorem that you can think of, anywhere in mathemat- 
ics? 

4. Why are there just dots in Figure 23.7, instead of curves? 

5. Suppose you're solving a system x' = Ax and find that the eigenvalues are 4.2 ± 3i. 
Which phase plane applies? Is the motion toward or away from the origin? How can you 
determine whether the motion is clockwise or counterclockwise just from the differential 
equation itself? 

6. What's rong with this? Figure 14.2 and Figure 24.1 don't match any of the 7 types 
listed in our classification theorem. So there must be a mistake in the theorem. 
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24 Nonlinear Systems in the Plane 



Today: Nonlinear systems in the plane, including things which cannot 
happen in a linear system. Critical points. Limit cycles. 

We now know about all the possible phase planes for systems of two linear 
equations. The next step is to see how this helps us analyze more realistic 
systems. There are two new ideas which enter into this study, and the first 
one is well illustrated by an example we looked at previously, the predator- 
prey system of Lecture 15. We will see that the phase portrait for the 
predator-prey system looks as though it has been made by selecting two 
of the linear phase planes from Lecture 23 and pasting them together in a 
certain way. In this sense, the linear systems prepare us pretty well for some 
of the nonlinear systems. Then we will show a second example containing 
entirely new behavior which can't happen at all in a linear system. So from 
this perspective it seems that perhaps the linear systems don't prepare us 
well enough at all for some of the nonlinear ones. 

We begin with the predator-prey system 



which we have discussed in lecture 15. The phase plane for this system is 
shown in the figure. 



x = x — xy 



y' = -y + xy 




-1 



-0.5 



0.5 



1.5 



Figure 24.1 The predator- prey system. Of course you have already noticed 
the part which looks like a center, and the part which looks like a saddle. 
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You may recognize some of its elements immediately if you have done the 
homework problems for Lecture 23 already. It shows a cycling between the 
sizes of the mouse and wolf populations. 

The points about which the action seems to be centered are known as critical 
points, and you may notice in Figure 24.1 that these are located at (0,0) and 
at (1,1). There is an easy way to find the critical points for any system, and 
that is to notice that they are just the constant solutions to the equations. 
In other words, all you have to do is set x' = and y' = and solve the 
resulting equations. 

Next we modify the system to include a little bit of emigration, thinking 
that perhaps the mice have discovered the inadvisability of living in this 
neighborhood. The system is now 

x = x — xy — .2 
y = -y + xy 

Observe the interesting effect this emigration has on the phase plane. What 
do you think it means, when that solution curve spirals outward for a while 
and then runs into one of the coordinate axes? It is not usually significant 
when a curve crosses an axis, but here we are talking about populations. 
When a population reaches zero, that is significant. The result of the emi- 
gration can therefore be extinction, or at least absence. 




Figure 24.2 The predator-prey system with emigration of the prey. Note 
the part which looks like a saddle and the part which looks like a spiral. 
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In both of these situations we have shown negative values of x and y. This is 
because the system exists independent of its interpretation as a population 
model, and we wanted to look at the whole phase plane. If your interest is 
in the population aspects though, you don't care about the negative values, 
and the conclusion is that as soon as a solution curve touches an axis, the 
model no longer applies. 

Our second main example is a classic known as the van der Pol equation. It 
is important for several reasons, partly historical and partly because it shows 
behavior which is impossible in a linear equation. The historical aspect is 
that this was originally a model of an electrical circuit, whose understanding 
was important in radio communications. The equation is 

x" + (-l + x 2 )x' + x = 

You should notice that it looks just like our second order equations from 
Lecture 10, except that the coefficient of x' is not a constant. So our char- 
acteristic equation method will do no good here. It converts to the system 

x' = y 

V = {l-x 2 )y-x 

which is not of the form z' = Az, so our eigenvector methods will also do 
no good. So we let pplane make a picture for us to study. 
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Figure 24.3 The van der Pol system contains a limit cycle. This is completely 
new, significant behavior not possible in a linear system. 
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Notice what is happening in this phase plane. There are solutions spiraling 
outward from the origin. That certainly can happen in a linear system as 
well. There are solutions coming inward also. But what is new is what 
happens between these. There seems to be a solution which is periodic and 
which the others approach. In the radio circuit what would that mean? It 
must mean that there is a voltage or current which cycles back and forth 
repeatedly. Now you recognize that this kind of thing happens when you 
have a center in a linear system too, but there is a difference here. Here there 
is only one periodic solution. This means that this radio is going to settle 
down near this one special solution no matter what the initial conditions are. 
That is very important. Such behavior deserves a new name; it is called a 
stable limit cycle. If you run this machine in your laboratory, the only thing 
you will observe is the limit cycle. 

Problems 

1. Run pplane on the system 

x — sm.(x + y) 

y = sin(xy) 

How many saddles and spirals can you find, visually? 

2. Consider the system 

/ 2 

r = r — r 
6' = 1 

The first one is the logistic equation. Look back at Lecture 2 to remind yourself that all 
positive solutions approach 1. What do the solutions to this system do then, if you take 
(r, 6) as polar coordinates in the plane? Convert the system to cartesian coordinates and 
run pplane on it. Did you find a limit cycle? 

3. Study a predator-prey system in which the predators are moving out. What happens? 
This can be used as a model of a fishing industry, in which the predators and prey are two 
fish species, and "move out" means get caught. Are the results reasonable or unexpected? 
If the plan is to harvest the sharks so that the tuna population will increase, what would 
your advice be? 

4. Run pplane on the system 

x — cos(x + y) 
y = cos(xy) 

How many limit cycles can you find, visually? 

5. Find all critical points for the system 

x = y(l~x 2 ) 
V = x + y 

Sketch the phase plane and try to figure out what type of linear behavior occurs near each 
of them. 
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6. Find all the critical points of the van der Pol system. 

7. Try to sketch a phase plane which contains two saddles. You are not asked for any 
formulas, but just to think about what such a thing could look like. Remember the 
uniqueness theorem, that solutions cannot run into each other. 

8. Try to sketch a phase plane containing two limit cycles. 

9. What's rong with this? Jane's Candy Factory used the system 

x = 10 - .25a; 

y = 6 + .25a: - .2y 

The critical points are the constant solutions, so we solve = 10 — .25a; getting x = 40, 
and = 6+ .25a — fly getting y = 32. 

25 Examples of Nonlinear Systems in 3 Dimen- 
sions 

Today: Nonlinear systems in space. Chaos. Why we can't predict the 
weather. 

We have seen lots of linear systems now, and we have also seen that these 
provide a background for some, but only some, things which can happen in 
nonlinear equations, in two dimensions. Today we move to three dimensions. 
We will see that nothing you learned about two dimensions can prepare you 
for three dimensions! There are aspects of systems of three equations which 
can be understood using the prior knowledge, but there are also things called 
"chaos" , which pretty well suggests how different and intractable they are. 
This lecture is only the briefest kind of introduction to this subject, but is 
here because it is important for you to know about these wonders. 

We begin with an example in which your knowledge of plane linear systems 
does help you understand a 3D system. Consider the system 

x = —x 
y = z 
z = -y 

These equations may be read as follows. The equation for x does not involve 
y or z, and vice versa. In fact the x equation is a simple equation for 
exponential decay. The y and z equations are a system corresponding to our 
favorite second order equation y" = —y whose solutions are sines and cosines. 
We are therefore dealing with a center in the (y, z) plane. Consequently we 
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can think of this system as having a circular motion in the (y, z) plane 
combined with an exponential decay in the x direction. So the solutions 
are curves in 3-space moving around circular cylinders centered along the 
x axis, and which approach the (y, z) plane as time goes by. Any solution 
already in that plane remains there, circling. 

In contrast with that system let us look at the Lorentz system 

x = 10(y — x) 

y' = — xz + 28x — y 
8 

z = xy z 

y 3 

This very interesting system was made by Edward Lorentz who was a meteo- 
rologist at MIT, as a very simplified model of circulation in the atmosphere. 
The system is nonlinear because of the xz and xy terms, and is derived 
from much more complicated equations used in fluid mechanics. You may 
be interested to know that the weather predictions we all hear on the radio 
and TV are to some extent derived from calculations based on observational 
data and the real fluid mechanics equations. Here is a phase portrait for the 
Lorentz system. 




Figure 25.1 The Lorentz system, (figure made using the Matlab commands 
shown in problem 2) This figure does not show a phase plane, but rather 
a projection of a phase portrait in 3 dimensions onto a plane. How do you 
suppose it occured to Lorentz, to think of the effect of a butterfly on the 
weather? 
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You may notice that there seem to be two important critical points about 
which the solutions circulate. One finds that the solution turns a few times 
around the left one, then moves over to the right one for a while, then back, 
and so on. The only problem is that it is not easy to predict how many 
times it will circle each side before moving to the other! Lorentz came upon 
this fact as follows. He was computing this on whatever kind of computers 
they had in 1963, and noticed the pattern of random-looking jumps from 
one side to the other. These were very interesting solutions, so he reran the 
computation from a point part way along, to study it in more detail. In the 
process, he typed in as initial conditions the numbers which had been printed 
out by his program. These numbers were slightly rounded by comparison 
with the precision to which they had been computed. The result made 
history: the qualitative features were still there, but the solution in detail 
was very different, in that after a few turns the numbers of cycles around 
the left and right sides came out different, just because of rounding off a few 
decimal places in the initial conditions. That means that if you try to make 
a picture like the one above, it will look qualitatively the same, but will 
be different in detail. It also means something about weather prediction. 
Lorentz asked, "Does the Flap of a Butterfly's Wings in Brazil Set Off a 
Tornado in Texas?" This is a restatement of the effect of rounding off, i.e. 
slightly changing, the initial conditions. Everybody has known for years 
that weather prediction is an art, that you really can't predict the weather 
very well. Part of what Lorentz did was to explain why you can't predict the 
weather: you can never know enough about the initial conditions no matter 
how many weather stations you build, when the equations of motion are so 
sensitive to initial conditions. 

We also have an enormous shift in the way the universe is conceived, as a 
result of studying this and other chaotic systems. On the one hand there is 
the uniqueness theorem which says that the future is determined uniquely by 
the initial conditions. This is the Newtonian approach. On the other hand 
we now know about chaotic systems in which yes, solutions are uniquely 
determined by their initial conditions, but they may be so sensitive to these 
conditions that as a practical matter we cannot use the predictions very far 
into the future. Isn't that interesting? It is a really major idea. 

Problems 

1. Check out the java applet de which is described in Lecture 7 and run some chaotic 
systems on it. Several are built in as examples. 

2. Solve a 3D system in Matlab. You can't use pplane of course, but you can set up a 
function file to compute your vector field, like 
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'/, file ef.m 



function vec=ef(t,w) 

vec(l)=10*(w(2)-w(l)) ; 

vec(2)=-w(l)*w(3)+28*w(l)-w(2) ; 

vec(3)=w(l)*w(2)-8*w(3)/3; 

'/, end file ef.m 

Then give commands like 

[t, w] = ode45(>ef ' ,0,20, [5 7 9]); 

plot(w(: ,1) ,w(: ,3)) 

This will compute the solution to the Lorentz system for < t < 20 with initial conditions 

x 

(5, 7, 9), and make a plot. The symbol w was used here for y , so what was plotted was 

z 

x versus z. Naturally you can play with the plot command to get different views. The 
hardest part is understanding why w(:, 1) means x. You have to realize that w is a list of 
three lists, before it makes any sense. Change these files to solve a system of your choice. 
Note that you are not restricted to just three dimensions either. 

3. Change the plot command in problem 2 to plot x versus t. 

4. Change the ef .m file of problem 2 to solve the linear system described in the text. You 
will have to play with the plot command to get a nice perspective view, for example you 
can plot x versus y + 2z or something like that. Does the picture fit the description given 
in the text? 

5. What's rong with this? According to this lecture, if you have a system of 3 or more 
variables you can get chaos. And, according to Lecture 15, if you have a system of 2 
spring-masses, you get 2 Newton's laws or a system of 4 first-order equations. Therefore 
2 spring-masses are always chaotic. 
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26 Boundary Value Problems 



Today: A change from initial conditions. Boundary values. We find that 
we don't know everything about y" — —y after all. 

You have by now learned a lot about differential equations, or to be more 
specific, about initial value problems for ordinary differential equations. You 
have also seen some partial differential equations. In most cases we have 
had initial conditions. At this time we are prepared to make a change, and 
consider a new kind of conditions called boundary conditions. These are 
interesting both for ODE, and also in connection with some problems in 
PDE. 

First you have to know what a boundary is. It is nearly the same thing in 
mathematics as on a map: the boundary of a cube consists of its six faces, 
the boundary of Puerto Rico is its shore line, the boundary of the interval 
[a, b] consists of the two points a and b. The concept of a boundary value 
problem is to require that some conditions hold at the boundary while a 
differential equation holds inside the set. Here is an example. 

Example: 

y = — y 

2/(0) = 
V (2tt) = 

We are asked to solve a very familiar differential equation, but under very 
unfamiliar conditions. The function y is supposed to be at and tt. The 
differential equation here has solutions y = A cos(i) + B sin(i). We apply the 
first boundary condition, giving y(0) = = A. So we must take A = 0. That 
was easy! Now apply the second boundary condition, giving y(27r) = = 
Bsin(27r). Well, it just so happens that the sin(2-7r) is zero. So the second 
boundary condition is fulfilled no matter what B is. Answer: y(t) — Bsin(t), 
B arbitrary. 

Notice how different this example was from our experience with initial con- 
ditions, that there were infinitely many solutions. Just the opposite thing 
can happen too: 

Example: 

y = — y 

2/(0) = 
2/(6) = 
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This seems close to the previous example since only 2n has been changed to 
6. This begins as before with A = 0. Then it just so happens that the sin(6) 
is not zero. So unlike the previous case, we have to set B = 0. Answer: 
y(t) = 0. There is no other possibility. 

What is a physical interpretation of these problems? We know that the 
equation y" = —y describes an oscillator, a rock hanging from a Slinky, 
vibrating up and down at frequency If you require the condition y(0) = 
0, that just means you want the rock at the origin at time 0. That alone is 
not a real problem, because you can start the stopwatch when it is there. 
Since the frequency of oscillation is the rock will be back at the origin 
after 2ir seconds. If you then require the condition y(27r) = then there is 
no problem because that happens automatically. But if you require instead 
that y(6) = 0, then you are asking for the impossible. The rock takes 
2ir = 6.283 . . . seconds to get back, period. The only exception is the 
special function 0, signifying that the rock never moved at all. In that case 
it will certainly be back in 6 seconds, since it is there already. That is why 
there is no solution in the second case except for the solution. A more 
down-to-earth example of a boundary value problem: "You can go out with 
your friends at 8:00 but you have to be home by 12:00." 

There is a tendency to use x as the independent variable when one is dis- 
cussing boundary value problems rather than t because usually in practice 
it is position rather than time, in which one is interested. We'll do that in 
the next example. There is also the concept of eigenvalue in connection with 
boundary value problems. 

Example: 

y" — —cy Finding c is part of the problem too. 
3/(0) = 

v00 = o 

The "eigenvalue" here is c. We will assume we are looking for positive values 
of c to keep things simple, and a homework problem will be to analyze the 
cases in which c might be or negative. To solve this equation y" = — cy we 
recognize that we are again looking at sines and cosines: y(x) = A cos( v / cx) + 
B sin( v / ca;). That is where our previous knowledge of differential equations 
comes in. If you have any doubts about these solutions, go no further! Go 
back study second order linear equations a little more. You can't build 
a house on a foundation of sand. All set now? The boundary condition 
3/(0) = = ^4008(0) = A determines A — as before. Then we have to 
deal with the other boundary condition 3/(71") = = B sin( v / C7r). Certainly if 
5 = this is satisfied, and we have found a solution y(x) = 0. But we have 
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another parameter c to play with, so let's not give up too easily. Question: 
Is it possible that 

sin(y / C7r) = 

for some values of c? If so, B can be anything and such c will be our 
eigenvalues. Well, what do you know about where sine is zero? Certainly 
sin(O) = but if c = then we are back at y(x) = which is uninteresting 
by now. Where else is the sine zero? The sine is zero at n, for example. 
That gives ^fcm = 7r, or c = 1. There is a good solution to our problem! 
y(x) — Bsin(x), B arbitrary, c = 1. But there are other places where the 
sine is zero. For example sin(27r) = 0. That gives */cn = 2n, or c = 4. 
So there is another solution y(x) = Bsin(2a:), B arbitrary, c = 4. Also 
sin(37r) = 0, giving a solution y(x) — Bsin(3:r), B arbitrary, c = 9. In 
general we have solutions 

y(x) = B sin(na;) B arbitrary 
c = n 2 n = 1, 2, 3, . . . 

So the eigenvalues for this problem are 1, 4, 9, 16, 25, ... 

26.1 Significance of Eigenvalues 

Until now you have known eigenvalues as something connected with matri- 
ces. Suddenly the same word is being used in a different way, although 
in both cases there is (some operation) (something) equal to (the eigen- 
value) (the same something). There is an abstraction here, that the equations 
Ax = Xx and y" = — cy have something in common. We need to suggest 
that there is much more to eigenvalues than this. They are very important. 
You have gotten along somehow until now without knowing about them, 
but they have been all around you! 

For example, the color of light from your flourescent lights as you read this 
is determined by the frequency of the light, which is mainly a difference of 
atomic energy levels for the atoms in the light bulb. These energy levels 
are eigenvalues of something, not of a matrix, but of Schrodinger's equation 
which was mentioned earlier. The music that you might be listening to 
as you read this comes in pitches which are eigenvalues of some equation 
describing the guitar strings' motion. The room you are sitting in may 
occasionally vibrate a little when somebody drops something heavy, and the 
frequency is an eigenvalue of something too. We will discuss just some of 
these things in this course, including the guitar strings but not the atoms. 
There is a drum vibration model in Lecture 34. 

Problems 
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1 . Find all the solutions to the c < cases of 



y — —cy 
3/(0) = " 
2/(tt) = 



2. Solve 

y = -cy 

2/(0) = 2/(5) 

Finding the possible values of c is part of the problem, but to keep things simple you may 
assume that c > 0. 

3. Find all the solutions to 

// 

y — —cy 

2/(0) = 

2/(3) = 



4. Find all the solutions to 



y — —cy 
2/(0) = " 
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27 The Conduction of Heat 



Today: A derivation of the heat equation from physical principles. Two, in 
fact. 

Our survey of ordinary differential equations is now complete, and we have 
done some work with partial differential equations too. As you know, they 
are also somewhat different from ordinary differential equations. For exam- 
ple, you could make up an ordinary differential equation almost at random 
and find that many of the methods we learned will apply to it, and the 
software will solve it approximately and draw pictures of the approximate 
solutions. You could do that, although the motivation for doing so might 
not be clear. In the case of partial differential equations, almost nobody 
would just make one up, and the reason is that the equations people are 
already working on, which tend to be about things in the real world, are 
hard enough to solve as it is. There are no cut-and-dried software packages 
for partial differential equations either, since to some extent each problem 
needs its own methods. What happens is that rather than solving a lot of 
equations as we did with ordinary differential equations, we study just a 
few partial differential equations with many different initial and boundary 
conditions. So partial differential equations are about real things and have 
names to reflect this fact. The one in this lecture is called the heat equation. 

The heat equation looks like 

ut = au xx 

It is an abbreviation for 

du d 2 u 
-(x,t) = a— 2 (x,t) 

The equation concerns the temperature u(x, t) of a metal bar along which 
heat is conducted. You can imagine one end of the bar in the fire and the 
other in the blacksmith's hand, to emphasize that the temperature may 
change with time and with position along the bar. The number a is a 
physical constant which we will explain while we derive the heat equation. 

Now let's read the heat equation carefully to see whether it makes any 
sense or not. Suppose you graph the temperature as a function of x, at a 
particular time. Maybe the graph is concave up as in Figure 27.1 on the 
left. What do you know from calculus about functions which are concave up? 
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Isn't the second derivative positive? The heat equation contains the second 
derivative, and says that the more positive it is, the bigger the time derivative 
will be. What do you know from calculus about the first derivative being 
positive? Doesn't it mean that the function is increasing? So we have to 
imagine what that means, and conclude that the temperature must be rising 
at any point of the bar where the temperature graph is concave up. Similarly 
for the right side of Figure 27.1. If the temperature graph is concave down, 
the temperature must be decreasing with time. Does that seem correct? 
There is no hope whatsoever of understanding partial differential equations 
unless you think about things such as that. 



Figure 27.1 Temperature versus position along a metal bar. The arrows 
show the time dependence, according to the heat equation. On the left the 
temperature graph is concave up and the bar is warming, and on the right 
it is cooling. The little slab will be used in the derivation below. 

As you will see from our discussion, it is not necessary to think of a metal 
bar. It is also possible to think of the conduction of heat in the x direction, 
where this axis passes through a door or wall, or any similar situation in 
which heat energy flows in one dimension only. If you want to talk about 
heat flow along the surface of a sheet of metal you will have to use the 
two-dimensional heat equation 



If you want to talk about heat flow throughout a solid block of metal or 
though still air, you use the three-dimensional heat equation 



Our plan is to derive the one-dimensional heat equation by showing what 
physical principles are behind it, and by using what you might call elemen- 
tary mathematics, i.e. calculus. Later we will show a different derivation 
for the three-dimensional heat equation using the exact same physical prin- 
ciples, but more sophisticated mathematics, the Divergence Theorem. In 
both cases we need to point out that we are assuming some familiarity with 




U t = a(u xx + Uyy) 



U t = a(U XX + Uyy + u zz ) 
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physics for these derivations. On the other hand it is not so much a specific 
knowledge of physics as a willingness to think about heat, that you really 
need for the next few paragraphs. Actually solving the equation later does 
not make so many demands as the derivation does, but thinking about heat 
will help a lot there too. We begin. 

Consider a little slab cut from the bar between coordinates x and x + Ax, 
and suppose the density of the bar is p mass per length. We are going to 
account for the heat energy contained within this little slab. We have 

Ax = length of the slab 
pAx = mass of the slab 

physical principle # 1: There is something called specific heat, c, which 
reflects the experimental fact that a pound of wood and a pound of steel 
at the same temperature do not contain the same amount of heat energy. 
According to this principle we have 

cupAx = the heat energy content of the slab 

physical principle # 2: Heat flows from hot to cold, and more specifically 
there is something called conductivity, k, which reflects the fact that copper 
conducts heat better than feathers do. According to this principle we have 

—ku x (x,t) = the rate heat enters the left side of the slab 
ku x (x + Ax, t) = the rate heat enters the right side 

These two principles give us two different expressions for the rate at which 
heat energy enters the slab. This is wonderful! Whenever you have two 
different expressions for the same thing, you are about to discover something 
important. Equate them: 

(cupAx)t = ku x (x + Ax, t) — ku x (x, t) 

Then 

ku x (x + Ax, t) — ku x (x, t) 

CpU t = T 

Ax 

Finally take the limit as Ax — ► to get 

cpu t = ku xx 

This gives the heat equation, and shows us that 

k 

a = — 
cp 
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Figure 27.2 A region R inside the main block of material, with some heat 
energy flowing between R and the rest of the block. Energy conducts into R 
through its boundary surface S. There is no heat source or sink, chemical 
reaction etc., inside R. You can place R anywhere within the material. 

A second derivation may be accomplished in three dimensions by using the 
divergence theorem. We imagine heat conduction in a solid block of some- 
thing. Fix any region R inside the solid and account for heat flow into R. 
The heat content is JJJ R cpu dV and its rate of change is JJJ R cput dV. The 
rate of heat flow through the boundary surface S is JJ S kVu ■ ft dA. By the 
divergence theorem this last integral is equal to j]J R V • (kVu) dV. Thus we 



for all R. Now if a function happens to integrate to zero it means nothing, 
except that here it is true for all R. This turns out to be enough to insure 
that the integrand is zero, and we get 



1. Suppose the same temperature distribution is present in two metal bars at a particular 
time, where the first bar has Ut = 3u xx and the second has ut = 4:U XX , due to different 
material. Which bar cools faster? 

2. A new kind of physics is discovered on Mars. Their heat equation says ut = au xx — u. 
Do Mars bars cool faster or slower than Earth bars? 

3. Amazingly enough, we have discovered that the heat equation used by the Klingons is 
different from either Earth or Mars. They have ut = au xx +u, Ut = a(u xx + u yy ) + u etc. 
Sometimes their pizza doesn't cool off at all. Why? 



have 




CpU t = k(U XX + Uyy + u zz ) 



Problems 
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4. There is a way of thinking about continuous functions which is a little beyond the level 
of these notes, but you might be interested to hear it: Suppose continuous function / is 
not zero at some particular point a, f(a) ^ 0. Maybe it is 0.266 there. Then you can take 
some other points x sufficiently close to a that f(x) is near 0.266, say f(x) > 0.200. This 
will hold for any point x in at least some small neighborhood of a. Knowing this, we take 
R in Figure 27.2 to be located right in that neighborhood. Then what can you say about 
MnfdV? 

This shows that if these integrals are zero for all R, then / must be indentically 0. 

5. What's rong with this? The heat equation can't be right because we already learned 
Newton's law of cooling in Lecture 8 and it doesn't say anything about the temperature 
depending on x. 

28 Solving the Heat Equation 

Today: The meaning of boundary values. Steady state solutions. Product 
solutions. 

There are many, many solutions to the heat equation, by comparison with 
any ordinary differential equation. For example all constants are solutions. 
This makes sense if you think about what the heat equation is about; a metal 
bar can certainly have a constant temperature. There is another special class 
of solutions which are also easy to find. These are the steady-state solutions, 
i.e. the ones which do not depend on the time. Everyday experience suggests 
that there are such cases; consider heat conduction through the wall of a 
refrigerator. The kitchen is at a constant temperature, the inside of the 
refrigerator is at a constant temperature, and there is a continual flow of heat 
into the refrigerator through the walls. To find such solutions analytically we 
assume that it is a function of x only. What happens to the heat equation 
then? We have ut = and u x = u'(x), so the heat equation becomes 
u" = 0. That is an ordinary differential equation, and a pretty easy one to 
solve. Integrating once we get u' = a. Integrating again gives u = ax + b. 
All straight lines are the graphs of steady-state temperature distributions? 
Check by substituting into the heat equation, and you will see that these 
indeed work. Note that these straight-line solutions are a generalization of 
the constant solutions. 

So far we have not said much about boundary conditions or initial con- 
ditions. It is typical of partial differential equations that there are many 
possible choices of conditions, and it takes a lot of work to decide what the 
reasonable ones are. Let's say for now that we are interested in heat con- 
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duction along a bar of length / which is located in < x < I. Then one 
set of conditions which specifies a solvable problem is as follows. You may 
specify the temperature at the ends of the bar, and the initial temperature 
along the bar. For example here is one such problem. 



u(0, t) = 40 a boundary condition at the left end 

u(5, t) = 60 boundary condition at right end; the length is 5 

u(x, 0) = 40 + Ax the initial temperature 

You may consider whether these conditions seem physically reasonable or 
not. The bar is going to have one end kept at 40 degrees, the other at 60 
degrees, and the initial temperature is known. Should this kind of informa- 
tion be adequate to determine the future temperatures in the bar? This is a 
hard type of question in general, but in this case the answer turns out to be 
yes. In fact we can solve this particular problem with one of our steady-state 
solutions. 

Practice: Find the steady-state solution to this problem. 
Next consider the problem 



Here we have a bar which begins uniformly at 40 degrees, and at time zero 
somebody presses ice cubes against both ends of the bar and holds them 
there forever. What should happen? In this case the initial condition does 
not match the boundary conditions, so the solution, if it exists, cannot be a 
steady-state solution. In this case, it is not possible to just quickly reason 
out whether the problem has any solution, or whether it might have more 
than one. Physically it does seem as though the bar should gradually cool 
off to zero, probably faster at the ends where the ice cubes are pressed. 

28.1 Insulation 

There is another popular boundary condition known as insulation. What 
does insulation do? If you put a hot pizza into one of those red insulated 



ut = 



u 



XX 



the heat equation 



u(5,t) 
u(x, 0) 
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containers that the delivery people use, the pizza still eventually cools off, it 
just takes longer, the better the insulation is. So the key to understanding 
insulation is that we model it by the assumption that no energy passes 
through it. From our derivation of the heat equation, we know that the 
rate of energy conduction is a multiple of u x . You might want to reread the 
derivation on this particular point. Bottom line: The boundary condition 
at an insulated end is u x = 0. 

Another way to think about insulation is like this: heat flows from hot to 
cold. So for heat to flow, there must be a temperature difference. On the 
other hand if u x = somewhere, then the graph of temperature has slope 
0, and microscopically there is no temperature difference for small position 
changes there. 

Now what is the following problem about? 



This is a bar of length /, insulated at the left end and held at 60 degrees at 
the right end. The initial temperature is 40 degrees all along its length. 

28.2 Product Solutions 

There is another class of solutions to the heat equation which can be found 
by a method we have used many times before. We try a solution of the form 

rt 

u = ce . 

If you have read Lecture 3 then you have seen a less specific assumption about u. 

This trick worked the first time on first order ordinary differential equations, 
and indeed the heat equation is first order with respect to time. The trick 
worked the second time with systems of linear equations. This time, we 
need to allow c to be a function of x. Substituting into the heat equation 
ut = au xx gives rce rt = ac"e rt . These will agree provided that ac" = rc. 
Whenever you try something like this you have to take a few steps along 
the path to see whether it is going to lead anywhere. In this case we have 
run across something that we know, so things look hopeful. What we know 
is how to solve ac" = rc. It is a second order linear ordinary differential 
equation for c, and it is one of the easy kind with no d term. To simplify it 



u x (0,t) 
u(l,t) 
u(x, 0) 



Uxx 



60 
40 
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a little let's write r = —aw 2 , then what we need is c" = —w 2 c. The solutions 
are c = Asm(wx) + B cos(tra). Thus we have found many solutions to the 
heat equation of the form u = (As'm(wx) + B cos(wx))e~ aw *. We call these 
"product solutions" because they have the form of a function of x times a 
function of t. 

Example: Determine the values of w such that there is a product solution 
to 

u t — au xx 
having frozen boundary conditions 

u(0,t) = u(ir,t) = 

We compute u(0, t) = Be~ aw ' in the product solution above. This will be 
zero for all t only if B = 0. So far we have u = Asin(u>:r)e -£ " 1 ' Then we 
compute u{ix,t) = A sin(w7r)e~ al " J . This will be zero for all t if we choose 
u>7T to be any of the zeros of the sine function. These make w =1, 2, 3,. ...We 
have found a list of solutions: 

— an 2 t 

u„(x,t) = A n sm(nx)e an , n = 1,2,3,.... 

Problems 

1. Describe what the following problem is about. 

Ut — Uxx 

u x (0,t) = 
u x {l,t) = 
w(a;,0) = 400 

Determine whether there is a steady-state solution to this problem. 

2. Describe what this problem is about. 

Ut — u xx 
u(0,t) = 
u{l,t) = 
u(x,0) = 400 

Determine whether there is a steady-state solution to this problem. 

3. Find a product solution to 

u t 

u(Q,t) 
u(tx, t) 
u(x, 0) 

4. Find a solution to 

Ut = 5u xx 



— u xx 

= 

= 

= 3sin(a:) 
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u(0,t) = 
u(n,t) = 
u(x,0) = 






2sin(:r) 



5. Find a product solution to 



Ut = 5it ; 



'XX 



u x {0,t) = 
u x (n,t) = 
u(x, 0) = 3 



6. Determine whether there is a product solution to 



u t = 



u x (0,t) = 
u x (n,t) = 
u(x,0) — 





3 



7. Suppose that u\ and 112 are solutions to 



u(0,t) 
u(tv, t) 








Note that no initial condition has been specified. Set s = u\ + U2- Is it true that s is also 
a solution to this system? 

8. Find a solution to 



which is a sum of two product solutions. If you need help, look ahead to the next lecture, 
for which this problem is intended to be a preview. 

9. Plot the graphs of the product functions u„(x, t) = sin(n:r)e~ an * versus x for several 
values of t, for the cases n = 1 and n = 2. Describe what is going on in your graphs, in 
terms of temperature in a metal bar. 

10. In problem 7, suppose u(x,t) = e~ at y(x), i.e. an exponential function of t times a 
function of x. What boundary value problem, of the type considered in Lecture 26 must 
y solve? 



u{0,t) 
u(tv, t) 
u(x, 0) 




29 More Heat Solutions 



Today: More heat equation solutions. Superposition. 
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Today we are going to analyze a heat conduction problem which happens 
to be solvable using some of the product solutions found last time. Our 
emphasis is not on formulas, but on understanding what the solutions mean, 
and thinking about the extent to which this heat equation corresponds with 
our daily experience of temperature. 

Consider the problem in which a metal bar with frozen ends has temperature 
initially as in Figure 29.1. We want to find out how this temperature changes 
as time goes by. What do you think will happen? We know from daily 
experience that heat flows from hot to cold, and in fact a quantitative version 
of that statement was a major part of our derivation of the heat equation. 
So probably the hot area to the left will supply energy to the cooler parts. 



1.2 




Figure 29.1 Our initial conditions. Note the spot at which the bar is hottest. 
What do you think will happen to the hot spot as time goes by, i.e. will it 
cool off, will it move to the left, to the right, etc.? (figure made by the matlab 
commands x = :. 031415 : 3 . 1415 ; plot(x, sin(x)+(l/3)*sin(2*x)) 

Our problem is as follows. 

u(0,t) = 
u(ir,t) = 

u(x, 0) = sin(x) + - sin(2x) 

Before attacking the problem analytically let's just read it carefully to see 
what it is about. That may be an obvious suggestion, but it is also a way 
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to prevent ourselves from doing something idiotic. We see from the first 
line that we are solving the case of the heat equation in which the physical 
constant a has been set to 1 for convenience. We can't see what the domain 
is until we read the boundary conditions. The second two lines say that we 
are dealing with the ice cube case, where the bar ends are at and it on the 
x axis and are kept frozen. 



Figure 29.2 Ice cubes at the ends of the bar. You might be worried about 
heat escaping from all around the sides of the bar. If so, think about heat 
conduction through a wall rather than along a bar, or else think of the bar 
surrounded by perfect insulation so that the x direction is the only one of 
any concern. 

Finally the fourth line of the problem says that the initial temperature is 
the sum of two terms. Graphed seperately, they look like Figure 29.3. 




Figure 29.3 The initial conditions are a sum of these two terms, (figure made 

byx = 0: .031415:3.1415; plot (x , sin(x) ) ; hold on; plot (x, (l/3)*sin(2*x) ) ; 

hold off) 

If you did problem 7 of Lecture 27, you may suspect why we have included 
Figure 29.3. Solutions to the heat equation may be added to produce new 
solutions. If you didn't do that problem yet, go do it now. In mathematics 
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this is called linearity, as we have emphasized several times in these notes. 
In physics it is sometimes called "superposition" to crystallize the idea that 
two or more things are going on at the same time and place. 

Now we are ready to solve this problem. There are product solutions to the 
heat equation of the form u = As'm(wx)e~ aw 1 for all values of w. Perhaps 
the ones we need are the cases of the sine function which appear in our 
initial conditions. Let's try that. Our candidates are 

ui(x,t) = ^4isin(x)e~* 

and 

u 2 (x,t) = A 2 sin(2x)e -4 * 

Now check to be sure that these solve the heat equation and the boundary 
conditions. Do that now. Does it work? Ok. Notice that it doesn't matter 
what the coefficients Aj are yet, and that the sum 

u(x, t) = A\ sin(x)e~* + A 2 sin(2x)e" 4 * 

is also a solution to the heat equation with these frozen boundary conditions. 
The next step is to attack the initial conditions. We need to compare 

u(x, 0) = A-y sin(x)e° + A 2 sin(2x)e° 

from our tentative solution, with the desired initial condition 

u(x, 0) = sin(x) + ^ sin(2x) 

Now the question is, can this be made to work? Well, yes, it is really not 
hard to see at this point, since we did all the hard work already. You just 
take A\ = 1 and A 2 = ^. This makes the initial conditions work, and we 
have checked everything else. So we have a solution to our problem. It is 

u(x, t) = sin(x)e~* + - sin(2x)e _4t 
o 

The next step is to understand what this solution means and how it behaves 
as time goes by. The best way to do this is for you to sketch the graphs of 
the solution for various times, thinking about what happens to each of the 
two terms in the solution. If you do that yourself, and it is a good idea to 
do so, you should get something like Figure 29.4 for the individual terms, 
and Figure 29.5 for the solution. 
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Figure 29.4 The two terms in our solution, plotted for t — and t = .3 Notice 
how much each of these has moved, in particular that the term involving 
sin(2a;) has moved a lot more than the sin(x) term. Why is that? 



1.2 




0.5 1 1.5 2 2.5 3 3.5 

Figure 29.5 The solution temperature at times t = 0, .3, .6, .9, 1.2. Now you 
see that the hot spot moves toward the center. Why? 

If you study the solution formula and the pictures you will see how significant 
the n 2 is in the time exponential part of the product solution. 

Problems 

1. Think about what should happen to a metal bar which has an initial temperature 
distribution consisting of alternating hot, cold, hot, etc., in a lot of narrow bands along 
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the bar. Think of a bar which is red-hot in six places and cooler in between. Do you think 
it will take very long for the temperature to change? Solve 



lit — V>xx 

u(0,t) = 

u(n,t) = 

u(x,0) — sin(12a;) 



and sketch the solution for t = 0, .3, .6. 
2. Solve 

— Uxx 

u(0,t) = 



u(n,t) = 

u(x,0) = sin(a;) + .5sin(3a;) + .25sin(5a;) 



Sketch the initial condition. 

3. Solve 

u x (0,t) = 
u x (n,t) = 

u(x, 0) = cos(a:) + .5cos(3a;) + .25cos(5a:) 
Note that the initial temperature satisfies the insulated boundary conditions. 

4. Run the Heat Equation ID applet mentioned in Lecture 7. Try various initial and 
boundary conditions to develop some feeling for what this equation is about. 
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30 The Wave Equation 



Today: Traveling Waves 

We have now solved several heat conduction problems and are about to turn 
to another partial differential equation known as the wave equation. You 
can figure out what it is about. Notice that as we were working on the heat 
equation we solved a fairly large number of problems distinguished from 
one another by boundary and initial conditions, while the equation itself 
never changed. This is typical of partial differential equations, that there 
are relatively few of them which are considered to be important, and the 
conditions can be changed to permit these few important equations to apply 
to many different settings. For the heat equation we only dealt with heat 
flow on a finite interval < x < I. Our wave equation can apply to many 
things also, such as the vibrations of a guitar string, but we are mainly 
going to talk about waves on an infinitely long string, the whole x axis. The 
reason for this choice is that guitar strings can be solved by methods similar 
to those we used for the heat equation, and we prefer to introduce new ideas 
here. 

The wave equation is 

u tt = c 2 u xx 

U tt = C 2 (U XX + Uyy) 

U tt = C 2 (U XX + Uyy + U ZZ ) 

in one, two, or three dimensions, respectively. Note that it is second order 
with respect to the time, like Newton's law. In fact it is Newton's law in 
disguise, sometimes. We know from Lecture 10 that second order equations 
can oscillate, so that at least is reassuring. In one dimension u may represent 
the vertical position of a vibrating string. 



Figure 30.1 The one dimensional wave equation describes vibrations of a 
string, or less accurately, water waves. 
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In two dimensions or three, the wave equation can represent vibrations of 
a drum head or of sound in the air, and other things like electromagnetic 
waves. It is a sound and light show all by itself. If you have read Lecture 34 
on drums, you will have seen the equation u t t = (ru r ) r . That is one version 
of the two dimensional wave equation, written in polar coordinates for the 
case where there is circular symmetry. 

Let's read the one dimensional equation again. If u is position, what is u t t? 
Think about it. Did you think about it yet? It is acceleration. Also when 
you draw the graph of u as a function of x as in Figure 30.1, what does u xx 
represent? Just as for the heat equation, it is the curvature of the graph. 
The wave equation says that when the graph is concave up, the acceleration 
must be positive. Does that seem to fit with reality? Think about swinging 
a heavy jumprope. Better still, get a rope and try it. See whether you think 
the acceleration direction matches the curvature in this way. 

The wave equation for a string may be derived as follows. We apply Newton's 
law to a bit of the string, accounting for vertical forces. 



Figure 30.2 The tension in the string pulls at different angles on the two 
ends of any little piece of it, depending on the curvature. 

The tension T in the string has a vertical component. On the left it is 
approximately Tu x (x, t) since the slope u x is the tangent of the angle, which 
is near the sine of the angle if the angle is small. On the right it is Tu x (x + 
Ax,t). The difference of these must be the mass times the acceleration, so 



where p is the mass per length of the string. Dividing by Ax and taking the 
limit as Ax approaches gives 




Tu x (x + Ax, t) — Tu x (x, t) = pAxuu 



T 

—u xx = u tt 
P 
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which is the wave equation with c 




Now let's try to solve the wave equation, using our usual method. We try 
u = f(x)e rt in uu = c 2 u xx . It becomes r 2 f(x)e rt = c 2 f"(x)e rt . Cancelling 
the time exponentials and setting r = ac gives /" = o 2 /. Wow! Another 
easy second order equation has popped out. Again we find that we must 
remember those second order equations. Go back and review if you don't 
remember. We find solutions f(x) = e ax , f{x) = e~ ax , and linear combi- 
nations. For example if a = 1 we have found the solutions u = e z e c * and 
u = e~ x e ct . To understand what these are, we graphed the second one in 
Figure 30.3. 




-2-1 01 2345678 



Figure 30.3 The wave u{x,t) — e~^ x ~ ' \ graphed for various times. Note 
that to keep the expression x — ct constant as t increases you have to increase 
x also, so this wave is moving to the right. It could make a surfer cry. 

So, we have found product solutions to the wave equation of the form 

u = e a(x+ct) and u = e< x ~^ 

for all values of a, even complex. What does this mean? It certainly means 
that there are an awful lot of solutions which are functions of x ± ct. Now 
let's get crazy and generalize. Suppose you have a function 

u(x, t) = g(x — ct) + h(x + ct) 
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not necessarily built from exponentials at all. If you have read Lecture 3 
you have seen waves traveling in one direction there. These go left or right. 
Taking derivatives by the chain rule we find 

ut = —cg'(x — ct) + ch'(x + ct) 

u a = c g + c h 
u x = g'(x — ct) + h'(x + ct) 
u xx = g" + h" 

Woah! Look what happened there, we got utt = c 2 u xx without assuming 
anything about g and h except that they be differentiable. Doesn't that 
mean you can have waves of nearly any shape moving left and right? Yes it 
does. They are called traveling waves. 

Problems 

1. Show that cos(a; — 2t) and l/(x + 2t) 2 are solutions to u u = c 2 u xx for some c. Find c. 
Graph these waves for t = 0, 1, and 2. 

2. The function sin(a: + 3i) + sin(a; — 3i) consists of 2 traveling waves and solves 



Utt 


— QUxx 


u(0,t) 


= 


u(-K, t) 


= 


u(x, 0) 


= 2 sin (x 


u t (x,0) 


= 



Sketch the two traveling waves, both separately and together, and try to see why the sum 
has the boundary condition at x = 0, even though neither of the traveling waves does 
so separately. 

3. Use the addition formula for the sine to show that the function in problem 2 is equal 
to 2sin(a;) cos(3t). It therefore represents a vibrating string having frequency J^. Then 
solve 



Utt 


— QUxx 


u(0,t) 


= 


u(tt, t) 


= 


u(x, 0) 


= 2sin(2a;) 


u t (x,0) 


= 



Sketch the motion for this case, and convince yourself that this vibration produces a sound 
one octave up from that of problem 2, i.e. has twice the frequency. 

4. Continuing in the tradition of problems 2 and 3, now use initial condition u(x, 0) = 
2sin(3a;). The note you get is the fifth above the octave, musically. 

5. This problem is about using the wave equation to model an echo. Consider a function 
of the form u(x, t) — f(x + ct) + f(—(x — ct)) where that is the same / traveling in both 
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directions. Show that u satisfies the condition u x (0, t) — 0. Think of sound waves moving 
in < x and a wall at the point x = 0. Sketch a graph of u for several times using a 
function / of your choice, and try to convince yourself that what you have sketched is a 
representation of waves bouncing off a wall. 

6. Repeat problem 5 for u(x, t) = f(x + ct) — f(—(x — ct)). This time you should get the 
boundary condition it(0, t) = 0, which involves a different kind of bouncing at the wall. 
Try to see it in your sketch. Which of problems 5 or 6 is more like a string vibration, with 
the string tied down at the origin? Which one is more like water waves at the edge of a 
swimming pool? 

7. Try the Wave Equation ID and 2D applets, and explore for a while to get some feeling 
for how the solutions behave. Try various boundary and initial conditions. 

8. What's rang with this? The function u(x,t) — cos(i) sin(x) couldn't be a solution to 
the wave equation because it isn't a traveling wave. 





Figure 30.4 A wave traveling initially to the left is shown at six times, t 
increasing from top to bottom in the picture. It is reflected downward and 
reversed left-right when it bounces off a wall at the left side. Does this fit 
Problem 5 or Problem 6 better? 
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31 Beams and Columns 

Today: A description of some differential equations for beams and columns. 



Figure 31.1 A beam with some exaggerated deflection. 

A beam is a bar loaded transverse to the axis, while a column is loaded along 
the axis. Of course the orientation doesn't matter even though we think of 
a column as being vertical, 




'//. 



Figure 31.2 Left: A segment of the beam, with shearing force V and moment 
M acting. R is the radius of curvature. 

Right: A column buckled by an applied force. 



For either case suppose x axis left to right gives position along the bar, and 
we have a segment of the bar as shown with shear force V [Newtons] and 
bending moment M [Newton-meters] as functions of x. Also F [N/m] is 
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an applied loading. Let y(x) be the downward displacement of the beam 
measured to the center line say, and R the radius of curvature at each point, 
also measured to the centerline. In mechanics class you will use 'centroid' 
which is more correct than centerline. 

Geometric fact, not included in most calculus courses: y" = j^. They do 
tell you about concave up and down (up is shown for the segment!) but this 
quantifies it. 

Practice: Write a function whose graph is the top part of a circle, and 
verify the second derivative at the topmost point. Does it matter if you add 
some slope mx + b to your function? 

The hardest part of the analysis is to see that y" is also proportional to M. 
Here is how that works. 

For the bending shown, the top portion of the beam must stretch, and the 
bottom compress. Try it with an eraser. Mechanical fact: steel has the prop- 
erty that for small stretching and compressing, the stress a = (force/area) 
and the strain e =(length change) /(original length) are related by Young's 
modulus E = 210 x 10 9 Newtons/ (square meter) as 

a = Ee 

For aluminum and granite E is about 1/3 as much. 

Next geometric fact: let y\ measure the distance from the centerline of the 
beam, away from the center of curvature. So y\ = at the centerline, 
yi = +2 at the top of the beam if the whole size of the beam is 4 units 
top to bottom, etc. Then the geometric fact is that the length at level y\ 
changes by the factor (1 + ^-). Exercise: check that out. Thus e = ^. 

So, the bending moment at a particular cross section must cause the length 
changes, and we have stress 

pulling or pushing normal to the cross section. Then the moment is the sum 
of force times moment arm, 

M = If a( Vl ) yi dA = E (I y\ dA(h = Eld) = Ely" 

J J cross— section J J cross— section 

This concludes the hardest part. 



128 



For example a q x q square cross section beam has moment of inertia I = 
f^J^yf dy ldy2 = qy 12. 

The easier parts are to see that M' is proportional to the shear V, and that 
V is proportional to F, the downward load [N/m] applied to the beam. 

Practice: For M 1 sum moments about the left end point of the segment, 
giving -M(x) + M(x + Ax) + (Ax)V(x + Ax) = 0. For V' sum downward 
forces on the segment, giving — V(x) + V(x + Ax) + F(x)(Ax) = 0. Then 
divide by Ax and take limits as Ax — > 0. 

Putting all that together we summarize beams: (sign conventions may differ 
in various texts) 

Ely" = M, Ely'" = -V, Ely"" = F 

These are differential equations, and in addition they tell us how to interpret 
various boundary conditions. For example a cantilever beam embedded 
solidly into concrete at one end probably has y = y' = at the concrete 
end, and y" = y'" = at the free end. 

Two more applications: 

For a column which is deflected into a bent shape y(x) by an axial force P, 
you get a bending moment y(x)P at each cross section, and taking account 
of signs the DE becomes 

Ely" = -Py 

Finally we mention vibrating beams: Suppose y depends on x and t and 
that there is no applied load F. If the mass of the beam is p [kg/m] then 
you find a PDE 

EIy xxxx = -py u 

for the vibrations. 

Problems 

1. Find a type of boundary value problem in Lecture 26 which could be applied to the 
column. Take -J^ = 1 if that helps. What is the smallest eigenvalue at which the column 
buckles? 

2. Try to find some beam vibrations of the form 

y(x,t) = acos(fti) sin(er) 

That does not mean 'derive this formula from scratch'. It means plug it into the beam 
equation 

PVtt = -EIy xxxx 
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and see what must be true of the numbers a, b, and c if this is to be a solution. 

But: it is alright to try a separation of variables too, as one way of deriving it. We are 
allowed to try more than one thing. 

3. Find a list of beam vibrations from Problem 2, which have the boundary values 

2/(0, t) = 0, y(n,t) = 



4. Use your list in Problem 3 to find a solution to 

ytt — yxxxx 

y(0, t) = 

V(«,t) = 

y{x,0) = sin(a;) - sin(3:r) 



5. The minus sign in the beam equation py t t = —EIy X xxx looks strange if you have been 
solving the heat equation recently. Suppose that at some time the shape of the beam is 
y = x 4 . Do you think the acceleration ought to be toward the +y or the — y direction? 

6. We know that the first overtone of a string vibrates twice as fast as the fundmental. 
What about a beam? Verify that it goes 4 times as fast when vibrating in two equal parts. 
This is why guitars and xylophones don't sound the same. 
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32 Power Series 



Today: How to make ugly functions beautiful. 



32.1 Review 



Polynomials are nice aren't they? 



f(x) = l, g(x) = 2-x 2 , 
h(x) = x 6 - 3.2x 3 + 4, y(x) = (200a; - 1) 
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They are the simplest functions you can build by repeatedly adding and 
multiplying numbers and one variable. People can use these as soon as they 
are have some idea of a symbol 'x' being used to represent a number. These 
are beautiful functions, and we know a lot about them. 

But we also know that there is not enough variety in those to do science. We 
need fractional powers, exponentials, logarithms, and trigonometric func- 
tions at least. With all those, we have a list of functions just long enough 
to include most of the buttons on a calculator. In other words, those are 
the functions which are so commonly used by every working scientist that 
it is worth the time of some manufacturing company to create a profitable 
product, the calculator, which includes them. They are not polynomials. 

Other functions are needed too, many of them actually defined by differential 
equations, that don't have their own button yet, so to speak. Maybe by 
the year 2035 the solution to some new equation will be so important that 
everybody has a calculator for it, built into their hat or something. 

All these functions wish they were beautiful too, like the polynomials. It 
turns out that some of them are nearly polynomials after all, though they 
may not know it. For example, it is sometimes possible to give pretty good 
approximations by polynomials. Here is one. 



cos(x) = 1 - -x 2 + —x 4 + E 



where \E\ < 




That is a nice result. For example, 



cos(.2) = 1 - .02 + 



.0016 



1 



= 0.98007 to four decimals. 



± 



24 



46080 
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You might remember doing such estimates in your calculus course under the 
topic: Taylor's Theorem. Sometimes it seems disturbing that you don't ever 
find out exactly what E is, because we are more accustomed to equations 
than to inequalities. But if you know E exactly, don't you know the cosine 
exactly too? And that is what we don't know. 

Functions that can be expressed as power series are known as analytic func- 
tions. Like the cosine, they can be approximated by polynomials not only 
of degree 4 but of any degree and to whatever accuracy you need. Here is 
what some power series look like: 

f(x) = c + c\x + c 2 x 2 + c 3 x 3 H 

y(x) = 2 + x - 1.7x 2 +4.3x 3 - 5x 6 

g(x) = co + ci(x-10)+c 2 (x-10) 2 +c 3 (x-10) 3 + --- 

e x = l + x + -x 2 + -x 3 + -x 4 + --- 

Lets look at both the form and the meaning of these expressions. 
The Form 

The first and third examples are generic ways to write a power series. The 
first is said to be expressed in powers of x, while the third is in powers 
of {x — 10). The coefficients ct represent various real or perhaps complex 
numbers. The 10 can be replaced by any other number, and the letter x can 
be replaced by any other letter, but no other change is allowed. For example, 
^fx does not occur because we only allow nonnegative whole number powers. 

The second example is a polynomial. That is a power series because it is in 
the form of the first example having most of the c& = 0, which is alright. 
But the other way around, a power series not considered to be a polynomial. 
Polynomials have a finite degree, the largest exponent which occurs. Also if 
you have an expression involving, say, (x 2 — 2) 3 , that is not a power series 
either, but it might be possible to rearrange it as one. For instance, in the 
identity 

( x 2 - 2) 3 = x e - 6x 4 + 12x 2 - 8 

both sides are polynomials, but only the right hand side is considered a 
power series. 

The fourth example is known as the power series for the exponential function, 
and we will say more about it later. For now just notice that whatever the 
meaning is, the equation seems to be giving us a beautiful almost-polynomial 
for e x . 
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Practice: Did you ever have a formula for e x before? Did your formula 
allow you to calculate, say e 0,573 without a calculator? 

So far we see what a power series looks like, but we haven't said a word 
about what you are wondering: what does the dot-dot-dot mean? How can 
we calculate anything that goes on forever? 

32.2 The Meaning of Convergence 

I'm going to explain the idea of convergence of power series by using one 
principal example. You can find a more traditional approach in your calculus 
book. 

Everything we know about power series comes from one example, known as 
the geometric series. It is the most important. You need to know everything 
about it, and how to recognize it even in disguise. 

It looks like this. 

l + x + x 2 + x 3 + x 4 + x 5 H = — — if \x\ < 1 

1 — X 

One disguised version of it occurs just below. As a first observation, note 
the restriction that |x| < 1. Most formulas you are familiar with have had 
few restrictions on the variable. But that needs to change now, and we'll 
have to think about the applicability of our formulas. For example in the 
geometric series, if x = 1 you would have complete nonsense on both sides 
of the equation, so that has to be excluded. Again if x = 2 for example the 
right side is defined but the left says add the numbers 1, 2, 4, 8, 16 etc and 
never stop adding. Obviously that doesn't work either. But lets try a small 
number, 0.1. Then it says 

1 + 0.1 + 0.01 + 0.001 + 0.0001 + 0.00001 + • • • = ^ 
Could that be true? It ought to be interpreted as 

1.1111111111 - • ■ = — 
9 

which seems to be correct. 

Practice: If you aren't sure about that, divide 10 by 9 and see what you 
get, both by hand and on a calculator. 
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We still haven't assigned any exact meaning to the dot-dot-dot expressions. 
Here we will first find out exactly what it means for the geometric series. 
We start with an algebraic identity. 

(l + x)(l-x) = l-x 2 

(1 + x + x 2 )(l - x) = 1 - x 3 
(1 + x + x 2 + x 3 )(l - x) = 1 - x 4 

Practice: Are those true? Is there any restriction on the values of x? 
In general we have 

(1 + x + x 2 + x 3 + • • • + x n ){l -x) = l- x n+1 
and so, if x / 1, 

\ — x 1l+l 

1 + x + x 2 + x 3 H Yx 11 = 

1 — x 

Now that is still a finite sum. Notice that if n is 1000, you have a very long 
expression on the left hand side, much too long to write out. But the right 
hand side is very short to write and pretty short to compute no matter what 
n is. What happens when we try to allow n — ► oo? On the left side, we 
don't have a clue because that is really what we are trying to define. After 
we make sense out of it, then it will be alright to write the left side as 

1 + x + x 2 + x 3 H 

but we don't know what it means yet. So look at the right side. What 
happens to 

1 - x n+1 
1-x 

as n — > oo? 

Practice: Try various numbers x, and convince yourself that there is no 
limit unless \x\ < 1, and then the limit is 
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So we understand the meaning of the geometric series as follows: If you stop 
after adding n + 1 terms you have 

1 + x + x 2 + x 3 H hx" 



1 — x 1 — X 
which is not exactly 



1 — x 

but the error is 

x n+i 



1 -x 

which approaches as n — > oo. 

Practice: What error does that give for the approximation 
1 + 0.02 + 0.0004 = 



1 - 0.02 

In general we make a definition, that the meaning of any series 

a + ai + a 2 + a 3 H 

is as follows, power series or not. Suppose there is a number L so that 
whenever you stop at the n-th stage you have some remainder R n 

ao + a\ + 02 + G13 + • • • + On = L — R n 

and R n approaches as n — ► oo. Then we write 

a + ai + a 2 + a 3 H = L 

and we say that the series converges to L. I hope you remember from 
calculus that ao + a\ + a 2 + 03 H h a n is called a partial sum of the series. 

In the next lecture we will see many uses of the geometric series. 
Problems 

1. Express the repeating decimal 

a = 0.306306306306 

as a (multiple of a) geometric series. Express a as a rational fraction a — m/n with 
integers m and n. 

2. More challenging: Consider the series 

_ J_ J_ J_ J_ 

S - + 2 2 + 32 + 42 + 52 + " ' 
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There is no number x so that the terms in s are each less than the terms of the geometric 
series 

= l + x + x 2 +x 3 -\ 

1 — X 

Try it, and you'll see why. But you can group the terms in an interesting way: 

• = 1 +(^ + ^) + (^ + - + ^) + (p + -+i^)+- 

so that every grouping is less than the corresponding geometric term with some carefully 
chosen number x. Try to figure out x. What can you conclude about the size of the 
number s? 

This series was a research problem somewhere around the 18th century. Eventually Euler 
found that it is a special case of a Fourier series, and he was able to work out the value 

2 

exactly, On page 167 you will work that out yourself. Even now you can see that 
there is something advanced about it, from the tt 2 . Nobody squares tv in everyday life. 

It is common in calculus books to condense this interesting story into the statement that 
'the ratio test fails for this series'. 
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33 More on Power Series 



Today: We use the geometric series as the spring from which a river of 
functions flows. 

We have learned that 

l + x + x 2 + x 3 H = (M < 1) 

1 — x 

and we have reviewed what the idea of convergence means. We have not 
been very quantitative yet about other series. And we have not justified 
why this is so important. 

So: why on earth would anybody want to replace the perfectly simple ex- 
pression 1/(1 — x) with a complicated limit, which doesn't even work for 
every value of x? 

I'm so glad you asked. Let me give you eight excellent reasons. 

1. We will quote a theorem below which says more about the convergence 
of power series, but for now lets just integrate our geometric series from 
to x. We get 

x + -x 2 + -x 3 + \x 4 H = - ln(l - x) if |x| < 1 

2 3 4 

Did you ever have a formula for the natural logarithm before? What was 
it? Did it allow you to calculate ln(0.9)? Lets try that: 

-ln(0.9) = 0.1 + - (0.01) + = (0.001) + • • • 

= 0.1 + 0.005 + 0.0003333- •• 

So ln(0.9) = -0.1053---. Nice, eh? Look Ma, no calculator! We have 
turned the ugly In into a beautiful almost-polynomial. 

2. Next make a change of variables in the geometric series to get 

l-u 2 + u A -u 6 = - — - — ^ (M<1) 

1 + u z 

and integrate from to x to get 

x-\x z + \x b = tan _1 (x) (|x| < 1) 
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Gee. We never had that before did we? 

3. As a limiting case in item 2, take x — > 1 to see 

111 7T 

1 -3 + 5"7 + --- = 4 

We have always been told what a strange number ir is, irrational, nonalge- 
braic, transcendental, infinite nonrepeating decimal expansion. People have 
written books about the hidden messages of the universe in the mysterious 
decimals of ir and other such nonsense. Here we see that somehow it is 
related to the odd integers. That is unexpected. 

4. Here is a standard textbook example: Determine whether the series 



OO j 

n=0 



converges. Lets write it out 



1 J_ , J_ 1 
+ 1 + 1-2 + l-2-3 + l-2-3-4 + "' 
The biggest single question is whether this could be finite, because of a 
fundamental fact about the real numbers. Write the n-th partial sum as 

1 1 1 1 1 

s n = 1 + 1 + 



2 2-3 n! 

Theorem: An increasing sequence of real numbers bounded 
above has a limit, i.e. if 

si < S2 < S3 < • • • < L < oo 

then the have a limit which is less than or equal to L. 

It is certain that our si < S2 < ■ ■ ■ • What about some bound L? Here is a 
way to check our sequence. Since 3, 4, 5, etc are all larger than 2 you have 
an inequality for each partial sum 

You see the geometric series lurking there in powers of |? So we can see 
that all of our partial sums 

Sn < 1 + T~ZA. = 3 
1 _ 2 
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That is enough to establish convergence (to something) by the theorem. It 
does not tell us exactly what the sum is, just that it is less than 3. Of course 
we are supposed to know it is e = 2.718 • • • but that comes from item 5 next. 
Also you can compare this approximation of e to the one on page 33. 

I hope you see that a knowledge of the geometric series gives information 
about many other series. In fact, if you have had a course in which the 
various "convergence tests" were discussed, it will be important for you 
to realize that most of those tests, such as the ratio test, are based on 
comparison with a geometric series just as in this example. 

5. The same kind of discussion as in item 4 reveals that the series 

f(x) = 1 + x + ^x 2 + ix 3 + ^x 4 + • • • 

converges for any value of x. Again it does not tell us that it converges to 
e x , as was claimed at the beginning of this section. We need a way to check 
that. We are going to use the uniqueness theorem for differential equations 
to do that. It is on page 22 if you need to look. Note first that the initial 
condition is 

/(0) = 1 

so it has that in common with e x at least. What characterizes e x l Isn't it 
the derivative? If we are allowed to calculate f'(x) from the series we get 

f'{x) = + 1 + ±2x + ^3x 2 + ^4x 3 + ■ • ■ 

= l + x + ^x 2 + -x 3 + --- = f(x) 

That is exactly what we expect of e x . We have the uniqueness theorem for 
solutions of differential equations, so it looks strongly as though f(x) = e x . 
In fact it is true. There is only one fault in the logic. We still need to know 
whether it is alright to differentiate the terms of a power series to calculate 
its derivative. This is similar to the problem we had in item 1 above, where 
we didn't know whether it was ok to integrate a series. Here is the theorem 
we need, from another course. 

Fundamental Theorem of Power Series For each power 
series 

oo 

^ ^ CnX 

n=0 
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there is a number R, called the radius of convergence. If —R < 
x < R, the series converges. If x < —R or R < x, the series di- 
verges. In fact the coefficients c n are allowed to be complex num- 
bers and the series converges for all complex numbers x inside a 
circle of radius R in the complex number plane, and diverges for 
every number outside that circle. 

Within the circle of convergence it is allowed to differentiate and 
integrate the series term by term, and the resulting series still 
has the same radius of convergence. 



The special cases R = and R = oo are used to announce that the series 
converges only for x = or converges for all numbers, respectively. The 
radius of convergence can usually be found by comparing the series with a 
geometric series, or by using any of the convergence tests, most of which are 
proved by comparison with a geometric series. 

6. Now that we have the e x series and the Theorem we get a lot of others 
almost free: 

cosh(x) = \{e x + eT x ) = 1 + ^x 2 + ^x 4 + ■■■ 

sinh(x) = -(e x — e~ x ) = cosh(x) = x + ^-x 3 + ^-x 5 + ■ ■ ■ 
2 ax 3! 5! 

7. And if we are comfortable with complex numbers we also get the trig 
functions: 

cos(x) = i( e « + e -») = 1 - ix 2 + ix 4 - • • • 

sin(x) = cos(x) = x — -x 3 + —,x b — • • • 
ax 3 5! 

8. How does a calculator work? Some of the functions it has buttons for 
can be calculated by series. Perhaps the buttons cause a partial sum of a 
series to be computed. Is that what it really does? Maybe you know, if you 
are in computer engineering. 

But what about functions that we don't know yet? The ones needed in 
science but not used enough to make them profitable as consumer products, 
as buttons on calculators? That is what the next section is about. 

Problems 
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1. What is the radius of convergence of 

1 + x + x 2 + x 3 + ■ • •? 

(Whatever you do, don't apply the ratio test to this series. The ratio test is proved based 
on what we already know about this series.) 

Does it converge when x = — |? To what? 

2. What does the power series theorem give for the radius of convergence of the derivative 
series 

l + 2x + 3x 2 + Ax 3 + ■■■ = (1 + x + x 2 + x 3 + x 4 + ■■■)'? 
Does it converge when x = — |? To what? 

3. What is the radius of convergence of 

1 + (x - 10) + (x- 10) 2 + (x - 10) 3 + ■ ■ ■? 
Does it converge when x — 8.3? 

4. Use your knowledge of the cosine function to figure out the value of 

_2 4 6 

IX TV 7T 

'-Y + Ji "61 + "' 
It really wouldn't pay to evaluate the first few terms of that, would it? 



34 Power Series used: A Drum Model 



Today: We set up a model for the vibrations of a drum head. 
The Drum 

Our model of a drumhead has radius 1, and it only vibrates in circular 
symmetry. That means you can only hit it in the center, so to speak. Using 
polar coordinates, there is dependence on r and t but not on 9. We write 

u(r,t) 

for the upward displacement of the drumhead out of its horizontal resting 
plane, and assume that u is small and that 

u(l,t) = 

since the edges of the head are pulled down against the rim of the drum. 

We will assume that there is a tension T [N/m] uniformly throughout the 
drumhead. Convince yourself that this is reasonable as follows. Imagine 
the tension forces on a small triangle of material located anywhere in the 
drumhead. 



141 



Suppose a triangle with vertices (0,0), (a, 0), and (0,6) has a force applied 
uniformly along each edge, perpendicular to the edge, in the plane of the 
triangle, with magnitude T times that edgelength. Show that the triangle is 
in equilibrium, i.e. that the sum of forces on it is 0, assuming T is a constant. 

Now write Newton's F = ma law for vertical forces on any annular region of 
the drumhead. We take the annulus to have radius running from r to r + dr. 
In working out the forces due to the tension T, it is best if you have already 
worked on the wave equation for string vibrations. The reason is that while 
the forces here are quite similar to those acting on the string, the geometry 
here is a little more involved because of the curved shape of the annulus. 




Figure 34.1 The annulus of drumhead material viewed from above. 

The length of the inner curved edge of the annulus is 2irr. There are forces 
directed inward all along that edge. We are interested only in the vertical 
resultant of those forces. That is T{2nr) times the sine of the small angle 
with the horizontal plane. The slope is u r (r,t). The slope u r is nearly the 
same as the sine of the angle for small displacements. See Problem 1 if you 
aren't sure about that. Similarly for the outer edge. 
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r 



r+dr 



Figure 34.2 A cut through the annular region, viewed on edge. The ver- 
tical displacement is greatly exaggerated. We're accounting for the vertical 
components of the tension forces. 

The net force vertically is then 

T(2vr)(r + dr)u r (r + dr, t) - T(2vrr)u r (r, t) 

The acceleration of the annulus is utt- Mr. Newton says that we also need 
its mass. The mass of the segment is its area times its density p [kg/m 2 ], or 
p(-7r(r + dr) 2 — 7rr 2 ) By Newton then 

T(2vr)(r + dr)u r {r + dr,t) - T(2irr)u r (r,t) = pir(2r dr + {dr) 2 )u tt 

Divide by dr and take dr — > getting 9 ^ T g^ = prutt- 
So our equation is ruu = ^(ru r ) r or equivalently 

T ( 1 

utt = — u rr + -u r 
P V r 

We will find some solutions to this having u(l,t) = at the rim of the 
drum. The solutions will be vibrations of our drum model, and hopefully 
will resemble vibrations of a real drum. 

Practice: Polar coordinates are awkward at the origin. What should the 
slope tt r (0, t) be at the center, so that the shape of the drumhead is smooth 
there? 

Separation of Variables 

Maybe our wave equation has some product solutions of the form u(r, t) = 
R{r)T{t). Since this is supposed to represent music, lets try 

u(r, t) = R(r) cos(u)t) 
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or sine rather than cosine. Set that into the PDE to find 
-J 1 cos(wt)R = ^ ^R" + ^R'^j cos(ujt) 

We need 

R » + 1 _ R i + ^E R = 

r T 

For any such function we can find, there will be drum vibrations R{r) cos{u:t) 
and R(r) sin(aji) 

In the next Lecture we plan to find a nice solution R. However, it will be 
easier if we don't have to deal with the parameters all the time. Let 
a new variable x = ar and choose a to clean up the R equation. Write 
R(r) = R(f) = f(x). We have R'(r) = af'(x) and R"(r) = a 2 f"(x). The 
DE for R becomes 

«V" + %/'+^/ = o. 

Of course we set a = uj^J~^. Our new DE becomes 

/" + -/' + / = 

x 

Problems 

1. Suppose a very acute triangle has vertices (0,0), (1,0), (l,m) so that the slope is m. 

(l,m) 



(0,0) f 1 - ) 



Convince yourself that m is the tangent of the acute angle exactly, and the sine of the acute 
angle is nearly equal to m. You could fiddle with the geometry, or try your calculator, or 
try the power series for the sine and cosine, or think of something else to try. 

That is why we feel it is justified to replace the sine of the angles by the u r values in our 
derivation. Of course, if it turns out that the model doesn't act like a real drum then we 
have to rethink this decision. 

2. If you are not clear about the limit we found, try this. First tell what this limit is: 

lhn f(x + h)-f(x) 



Then how about 



Finally the one we need 



lim 



h 

(x + h)g(x + h) — xg(x) , 



^ (r + dr)u r (r + dr,t) — ru r (r,t) 
dr^o dr 
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35 A new Function for the Drum model, Jq 



Today: We learn a new function which might not be on your calculator, 
and use it to describe vibrations of the drum head. 

So far we have derived a wave equation 

T ( 1 

Utt = — Urr + -Ui 
P V r 

for the displacement of our drumhead, and we know that we need to solve 
f" + ^f'+f = 0. Having such a function /, we will get many drum vibrations 

u(r,t) = f(u^r) cos(ujt) 

A new function, Jo 

Try a power series 

f(x) = fo + hx 2 + / 4 x 4 + • • • 

If you did the Practice item on page 143, you know that we don't want a 
term f\x. We are guessing that we don't need any of the odd degree terms. 
You can check later that we truly don't. To keep things as clean as possible 
lets also assume /o = 1. We can rescale things later as needed. We have 

/" + V + / = (2/2 + 4 • 3/ 4 x 2 + 6 • 5/ 6 x 4 + 8 • 7/ 8 x 6 + • • • ) 

+ (2/2 + 4/ 4 x 2 + Qfox 4 + 8f 8 x 6 + •••) + (l + f 2 x 2 + / 4 x 4 + fax 6 + • • • ) 

= (4/ 2 + 1) + (4 • 4/4 + f 2 )x 2 + (6 • 6/ 6 + / 4 )x 4 + (8 • 8/ 8 + k)x & + ■■■ 
For that to be identically 0, we need all the coefficients to be 0, so 

f - _! f - h ^ 1 f ~-f±- -±- f 1 

J 2 r>9 ' A?. n9.A'?.i a"). o5^2e5' J 8 



2 2' ^ 42 2 2 4 2 ' 17 D 6 2 2 2 4 2 6 2 ' J 2 2 4 2 6 2 8 2 ' 

etc. So 



x 2 X 4 X 6 X 8 



1 r,9 + o5^2 05/1^/25 "I" 



2 2 1 2242 2 2 4 2 6 2 2 2 4 2 6 2 8 2 
Isn't that a beautiful series? This function is sufficiently useful that it has 
a name, even though you won't find it on every calculator. It is called the 
Bessell function Jq. 
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5 10 15 20 25 30 35 40 

Figure 35.1 Graph of Bessell J (x) = 1 - §^ + ^ - H 

So /(x) = Jo(x) or, if we need a multiple of that, f{x) = foJo(x) for some 
constant /o- We need to know what this function is like, to see whether it 
gives us plausible drumhead vibrations 

u(r, t) = /oJo(wJ^r) cos(wi) 

The series tells us some of the properties that we need, and the differential 
equation 

4 + - j + j = o 

itself gives other properties. 

First the series: The series is reminiscent of that for the cosine function, so 
maybe Jo has many zeros, and oscillates, and is periodic. Well, 2 out of 
3 isn't bad: Jo has infinitely many zeros and oscillates somewhat like the 
cosine function, but is not periodic. In a homework problem you will use the 
series to find that Jo (4) is negative. Since we know Jo(0) = 1, there must 
be a number x\ between and 4 where Jo(xi) = 0. In fact the first root x\ 
is roughly 2.5 and the second root X2 is near 5.5, X3 is about 8.5, and there 
are infinitely many others. 

From the DE we see this: suppose you have an interval like [xi,X2] where 
Jo is at each end, and nonzero between. At any local min or max where 
Jq = 0, the DE then tells us that Jq = —Jo there. The graph can't be 
concave up at a local max, nor concave down at a minimum. So the graph 
can only have simple humps like the cosine does, no complicated zig-zags 
between the roots. 
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After analyzing Jo like this for some time, people eventually published tables 
of its values, and many years later someone built it into matlab and octave 
so that we could look at the graph in Figure 35.1. Isn't that nice? 



35.1 But what does the drum Sound like? 

We don't know oj yet. We don't know what sounds the drum can make. 

We have two boundary conditions. The first is that u r (0,t) = 0, so that 
the drum shape is not pointy in the center. This is automatic since it is a 
property of Jo- 

Practice: How do you know from the series 

J (x) = l-^ + ... 

that Jq(0) = 0? 

The second boundary condition is that at the rim we have 
u(l, t) = foM^\ ^ • 1) cos(ut) = 



To acheive this, we need 



T 



P 

w< / — = Xi 
w T 



or £2 or some root of Jo- Then the lowest sound we hear will correspond to 
the case 

Ft 

u ( r , t) = foJo{xir) cos(W —x\t) 

The ordinary cosine function has period 2ir, so this one will go through one 
cycle as 

ft 

-X\t 

P 

goes from to 2ir. One cycle in every M seconds. The frequency is 

4 / — Xl 

VP 

p fjr ^ or l ar S er ro °ts X2 etc you get shorter cycles, higher pitched sounds. 
Any linear combination of these various solutions, for different roots x n , and 
using both cosines and sines, is a solution to our wave equation. 

So our model predicts what sounds the drum will make. Good. 
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Practice: When you tighten up the drum tension T, do the frequencies go 
up or down? Also the drummer might change to a heavier material so that 
p increases. What does that do to the frequencies? 



Example: Take p = T for simplicity. We have a vibration 

ui(r,t) = Jo(xir) cos(a;it) 

where x\ is approximately 2.5. This is a solution to the boundary value 
problem 

Utt = U rr + -U r (0 < T < 1) 

r 

u(l,t) = 

U r (Q,t) = 

having initial position ui(r,0) = Jo(xir). Another solution is 

U2(r, t) — Jo(x2r) svn{x2t) 

where X2 is approximately 5.5. This has a different initial shape, and a 
different frequency, and a different initial velocity. 

Practice: What is the initial velocity of ui? 
Problems 

f . In this problem you can confirm part of what is shown in the graph of Jo, that Jo (4) < 
so there must be a root x\ somewhere in the interval (0, 4). Write the partial sum of Jo (a;) 
through sixth powers as sq (x) 

x 2 x 4 x 6 
s e(z) = 1 - 22 + ^242 _ 524262 

First check that S6(4) is negative. It saves a little work if you notice that the second and 
third terms cancel. 

2. In Problem 1 the tail of the series is so small that Jo(4) itself is also negative. That 
means, 

Jo (4) = s 6 (4) + E 

where E is not very big and the sum is negative. You can estimate E by observing that 
the tail is an alternating series, if you remember those from calculus. 

3. Make a sketch of the graphs of Jo(a;ir) and Jo(a;27") for < r < 1, so that you can 
see what some of the drum waves look like. You can figure these out from Figure 35.1 by 
rescaling. 

4. For a particular real drum, it might be feasible to measure p if you could weigh a sample 
of the drumhead material, but T could present more problems. Suppose you compare the 
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sound of the drum to a piano or guitar and decide that the fundamental tone is near the 
note A 2 = 110 [cycles/sec], the next to lowest string on a guitar. Figure out j for this 
drum. 

5. Remind yourself that the first three tones of a string have frequencies in the proportion 
1:2:3, i.e., there are solutions cos(t) sin(x), cos(2t) sin(2x), and cos(3t) sin(3a;) to the 
wave equation y tt — y xx . But what about a drum? Check the graph of J and find out 
approximately how the lowest three frequencies of a drum are related. This is why pianos 
and drums don't sound the same. 
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36 The Euler equation for Fluid Flow, and Acous- 
tic Waves 



Today: A model for compressible air flow, including sound. This example 
is a little more advanced than the rest of the lectures, so there are several 
Practice items here to guide your reading. I bet you can handle it. 

The purpose of this lecture is to derive three PDEs which describe compress- 
ible flow of air in one dimension, say in a pipe. You might imagine pumping 
air into a bicycle tire through a hose, for example. The main assumptions 
are: 

• There is no conduction of heat through the air. There is no friction 
of the air with the pipe, or any exchange of heat between the air and 
the pipe. 

• You know about Newton's law F = ma, and about the idea of conser- 
vation of mass. And of course you know calculus. 

Two of the three equations are known as the Euler equations for homentropic 
flow. You do not have to know what 'homentropic' means. Then having the 
Euler equations, we can further derive that certain approximate solutions 
ought to also give us solutions to the wave equation. These solutions are very 
familiar-they are sound waves. Thus we also get, almost free, a derivation 
of the wave equation. If you have already derived the wave equation for 
vibrations of a string, you will see that there are two very different situations 
giving rise to the same mathematics. 

The method used in this derivation is rather common in science. There are 
many assumptions behind every statement. So, to understand the deriva- 
tion, you have to ask yourself at every step, Why should that be true? What 
is being assumed? This careful questioning can lead to good insights. The 
biggest differences between this lecture and what you might find in a more 
advanced treatment, are that 1) we work in one dimension rather than three, 
and 2) we make some of that careful questioning explicit. But you can still 
find plenty of additional questions to ask. 
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36.1 The Euler equations 



The pipe is full of air. The air velocity is u(x, t) [m/sec], and density p(x, t) 
[kg/m 3 ]. For air, the density is typically around 1.2 near the surface of the 
earth, but the velocity could vary over a large range. We can think of x and 
u positive toward the right. 



w 



Figure 36.1 A portion of air moving in the pipe. The pipe has cross-sectional 
area A. The ends L and R don't have to move at the same speed because 
this is compressible flow. 

We focus on a portion W of the air. The letters L and R in the figure refer 
to the left and right ends of the mass W. W is mathematically just a moving 
interval. The left coordinate L is moving at velocity ul, and to be more 
precise, 

L'(t) = u(L(t),t) and ul is an abbreviation for this 
and similarly for the right end R at velocity ur. 

Besides velocity u(x,t) and density p(x,t), we must keep track of the air 
pressure p(x,t) [N/m 2 ]. It is very easy to mistake p and p typographically, 
but this won't happen to you because you are reading very slowly and think- 
ing hard, right? The pressure provides forces pi,A and PrA on the sides of 
W. The homentropic assumption is that the pressure is related to the den- 
sity by p = ap" 1 , where 7 = 1.4 and a is constant throughout the flow. Aside 
from this thermodynamic assumption, the rest of our derivation will consist 
of ordinary mechanics. 

Newton: The rate of change of momentum of the matter in W is equal to 
the sum of applied forces. 

R 

puA dx = —PrA + p^A 
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Practice: Is this really F = ma or is it more like j(mi)) = Fl When 
does the m in Newton's law have to be included inside the derivative? Of 
the quantities and expressions p, u, A, and J dx, which constitute v and 
which ml Why do we multiply pressure by area? Why is there a minus sign 
on prA? 

In the left side of Newton's law we need to remember that 

d r m r f{t) 

- h(x,t)dx= ht(x,t)dx + h(f(t),t)f'(t)-h(g(t),t)g'(t) 

dt J g (t) Jg(t) 



Practice: Problem 1. Derive that from the Fundamental Theorem of 
Calculus and the Chain Rule. 



So, if we cancel all the A's, Newton's law and the homentropic assumption 
give 

R 

(pu) t dx + {pu) R u R - {pu) L u L = -p R +pl = ~ap R + ap\ 

Now we may be able to study a limit of this for very short intervals. If it is 
possible to divide by R — L and take L — > R, then we get 

(pu) t + (pu 2 + ap r ) x = 

This is a partial differential equation which of course assumes the various 
functions are differentiable. We assume this. But it is important to know 
that the limit has been found not in agreement with all experiments. We 
won't go into 'shock waves' at all, but you should know that there are such 
possibilities. 

We notice that we have two unknowns u and p. Probably one equation isn't 
going to be enough, so we look for another equation. Turning to conservation 
of mass, we figure the mass in W to be pAdx. This doesn't change 
as the fluid moves because our set W is moving with the fluid. We have 
conservation of mass 

R 

pA dx = 

Doing the same transformations as we did for Newton gives 
f R 

J [p t A dx - {pu) R A + pu) L A = 
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and in the limit 



Pt + (pu) x = 



Practice: Problem 2. Use 

p t + {pu) x = 

in the Newton PDE to find 

tit + uu x = —a'yp' 1 ~ 2 p x 

These two equations are the results of our work, the one-dimensional home- 
ntropic Euler equations. 



Practice: Problem 3. We want to identify a physical meaning for the 
combination of terms u t + uu x which occurs in the equation of motion. We 
can show that it is the acceleration of the particle which is passing through 
the point x at time t. Suppose a particle of fluid has position given by a 
function x(t). Then we have two expressions for the velocity of this particle, 
x'(t), and u{x(t),t). These must be equal. Differentiate to show that the 
composition (u t +uu x )(x(t),t) is equal to the acceleration of the air particle. 



Practice: Can there be any acceleration at a point where ut = 0? 

36.2 Sound 

Here we start from the Euler equations, and imagine a small disturbance 
superimposed over an ambient stillness. 

Practice: Is it true that the ambient stillness u — 0, p = (constant) is a 
solution to the Euler equations of Problem 2? 

We then look for approximate solutions of the form u(x,t) = ev(x,t), 
p(x, t) = p\ + ew(x, t) where e is supposed to be a small number. The 
constant p\ could be taken to be a typical air density at sea level on earth. 

Practice: Problem 4. Substitute our assumed u and p into the Euler 
equations and ignore everything containing e 2 or higher powers. Show that 
you find for the e 1 terms, 

Vt = —ajp1~ 2 w x 
w t + piv x = 
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Combining the two PDEs in Problem 4, we have wu = —piv x t = a-7Pi w xx 
or 

w t t = c 2 w xx 

where c 2 = a^ypj 1 So the pressure disturbance satisfies the wave equation. 

Practice: What did we just now assume about the second derivatives? 
This subtlety is omitted routinely in scientific discussions. Are solutions 
to an approximate equation necessarily approximate solutions to the right 
equation? This point is also usually omitted. 

Air pressure at sea level on the earth is about 10 5 [N/m 2 ]. Using this 
and previous information, estimate the value of the coefficient c 2 occur- 
ing in our wave equation. Then, recognizing that traveling waves such as 
f(x — bt) + g(x + bt) satisfy a wave equation, what do we learn about the 
speed of sound? If lightning strikes at a distance of 3 football fields from 
you, how long before you hear it? 

Finally, we have to ask: Have we proved anything about the behavior of 
real air? What is the criterion for correctness in Physics? in Mathematics? 
Have these notes strictly conformed to either? 

Loose ends 

a) We decided to look for 2 equations for our 2 variables u and p. How reliable is the idea 
that one needs n equations in n unknowns before you can solve anything? Have you seen 
examples in linear algebra, Ax — b, where two equations in three variables might have no 
solution? where five equations in three variables might have a unique solution? Think 
of the eigenvalue problem Ax — Xx, which is nonlinear because we view both A and x as 
unknowns. Suppose A is 8 x 8, how many variables are there? equations? solutions A? 
solutions xl There is some truth in this idea, but few guarantees. 

b) You might be wondering whether we need to have another equation about conservation 
of energy. Ordinarily we would, but our homentropic assumption p = ap 1 with no heating 
and no friction, has the consequence that all energy changes are accounted for already 
by Newton's law. For a simpler example, a frictionless mass on a spring has Newton's 
law mx" = — kx. After multiplying by x' you can integrate once to get the conservation 
of energy ^rn(x') 2 + \kx 2 = constant. A similar thing happens here. If we were to 
allow heating of the air, then we would need another variable which could be taken to be 
temperature, and another equation. 

There is plenty of room for misunderstanding about 'heating': the word heating means 
transfer of heat energy, as through the walls of the pipe. The temperature, which is not 
the same as the heating, does in fact change even without heating, because air obeys the 
ideal gas law p — 287 pT [mks] . 

d) We showed that the density disturbance w satisfies the wave equation. You can show 
similarly that the velocity disturbance v satisfies the wave equation with the same coeffi- 
cient c 2 . What is the significance of the fact that the velocity and pressure disturbances 
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satisfy the same wave equation? Specifically, is it physically plausible that these two 
disturbances might travel at different speeds? 

Projects 

Project 1. There are at least two ways to understand our earlier statement about the mass 
integral, that the mass in W is pA dx. You probably know that you can think of dx 
as a bit of length, Ada: as a bit of volume, pAdx as a bit of mass, and the integral adds 
them. Another interesting approach is explained by P. Lax in his calculus book. Suppose 
we write S(f, I) for the mass contained in interval / = [a, b] when function / gives the 
density in I. S has two physically reasonable properties: 

1. If I is broken into nonoverlapping subintervals I = 7i U I2 = [a, c] U [c,b], then 
S(f,I) = S(f,h) + S(f,h) 

2. If there are numbers m and M such that m < f(x) < M then m(b — a) < 
S(f, [a,b])<M(b-a) 

Under these conditions, you can actually prove that S(f, [a, b]) = f(x) dx. The project 
is to look up whatever definitions you need and figure out why this works. The idea applies 
to many applications other than mass. 

Project 2. Our derivations tracked the momentum and mass of a moving portion of 
fluid. Some people prefer instead to track the momentum and mass inside a fixed interval, 
accounting for stuff going in and out. Do you think one of these approaches might be more 
correct physically than the other? The project is to redo the Euler equation derivation 
by considering a nonmoving interval N = [a, b] instead of the moving W. The mass 
conservation is easiest: The rate of change of the mass which is currently contained in N 
is 

Jt[ PAdX 

and this ought to be equal to the rate in at b plus the rate in at a: 

= -pbAllb + p a Au a 

Of course this leads to the same mass PDE as before. Try to redo the Newton law in this 
context. 



155 



37 Exact equations for Air and Steam 



Today: Some vector fields are gradients, and some are not. So some differ- 
ential equations are said to be exact. Here is the history. 

It was in the early days of steam engines, when people first found out that 
there was a new invention on which they could travel at 25 miles per hour. 
No human had ever gone nearly that fast except on a horse, or on ice skates. 
Can you imagine the thrill? 



Figure 37.1 How fast can it go? 

It was an outgrowth of the coal mining industry, of all things. People used 
coal to stay warm, and unfortunately for the miners, the mines tended to 
fill with water. A pump was made to fix this problem and it was driven 
by an engine which ran on, well, it ran on coal! But people being as they 
are, it wasn't long before somebody attached wheels to the engine and they 
started competing to see who could go fastest. 

At about this time people noticed that every new train went faster than 
the last one. The natural question was whether there was any limit to the 
speed. So M. Carnot studied this and found that he could keep track of the 
temperature and pressure of the steam, but that neither of those was equal to 
the energy of the moving train. Eventually it was worked out that the heat 
energy added to the steam by the fuel was indeed related to the temperature 
and pressure. They called the new rule the first law of thermodynamics. It 
looked something like this, although the numbers I'm using here are for air, 
not steam: 



is supposed to hold whenever a process occurs that makes a small change 
in the temperature T [Kelvin] and the specific volume V [ra 3 /kg] of the 
gas. Here, the pressure comes in, again for air, through the ideal gas law 
P = 287 pT, V = 1/p. A main point discovered: that expression is not a 
differential. This was so important that they even made a special symbol 
for the heat added: d Q which survives to this day in some books. 




n 



=Qco 



heat added = 717 dT + 287 T 



dV 



[Joule /kg] 



156 



Practice: Using simpler numbers and variables, check that 

7dx + 2-dy^df 

y 

for any function f(x, y). If it were df = f x dx + f y dy, then you would have 

x 

7 = f x and 2- = f y 
V 

See why that can't be true? What do you know about mixed partial deriva- 
tives, f xy and f yx ? 

This relates to the gradient because an alternate way to express that is: For 
every function / 

7f+2-jV V/(z,y) 

y 

Anyway, the big disappointment to the steam engine builders was that the 
energy added in the process was not a differential, which meant that you 
could not make a table of values for how fast you are going to go, based only 
on temperature and pressure. 

But, the big discovery was that if you divide the heat added by T, you can 
make a table of that. 



Practice: 

7 — + 2 — = d(something) 
x y 

and so for air you also have 

„-„dT „„„ dV heat added ,, . 
717 — + 287 — = = d(somethmg) 

The "something" is called entropy, S, and for air, 717 +287 ^ = dS. In 
Problem 1 you can figure out S from that. There are tables of the entropy 
of steam in the back of your thermodynamics book. 

So what did Carnot come up with? Is there a limit to how fast the train 
can go? Well, he thought of an idealized engine cycle where for part of the 
time you have dT = and for the other part, dS = 0. Since it is possible 
to relate the V changes to the mechanical work of the engine, that allows 
a computation to proceed. He worked out how fast that ideal train can go. 
You'll have to read about it in your thermo book. 

Problems 

1. Do the Practice items if you haven't yet. Work out the entropy of air as a function 
of T and V. 

2. Using P = 287 pT work out the entropy of air as a function of T and P. 
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38 The Laplace Equation 



Today: Two roads lead here. The heat conduction and drum vibration mod- 
els both morph into the Laplace equation under certain conditions. Many 
roads lead away. The one we take will be the Fourier road. 

The Laplace equation looks like this in two dimensions 



or in one u" = 0, or in three u xx + u yy + u zz = 0. In two dimensions using 
polar coordinates not at the origin it looks like 



Solutions are called "harmonic" functions. In every case it says that the 
divergence of the gradient of u is 0, div(Vu) = 0. It is abbreviated An = 
by mathematicians and V 2 u = by engineers. But what does it mean? We 
approach that from two directions, as follows. 

I. The Laplace equation can apply to a steady state temperature. We know 
from Lecture 27 that the temperature in a region is modeled by the heat 
equation ut = Au, which keeps track of the conduction of heat energy from 
hot places to cold places. If the temperature is not changing with time, 
that gives the Laplace equation. It doesn't mean there is no energy flow. A 
refrigerator can maintain a constant temperature but only by forcing out the 
energy that seeps in through the walls. We will consider the polar coordinate 
case in a disk with the temperature u specified on the boundary circle. 



Figure 38.1 A membrane stretched, or the graph of a steady temperature. 

II. Our drum discussion in Lecture 34 can be modified. In that case we had 
a dynamic situation u rr + \.u r = utt with circular symmetry. If you instead 
analyze a circular membrane stretched statically over a varying rim height, 




u rr + -u r + -^uee = 0. 
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you find for small angles that the same Laplace equation u rr + -u r + ^uqq = 
does the job for the displacement u. What is added is the 9 dependence 
including nonzero boundary conditions, and what is removed is the time 
dependence. 

So we are solving the boundary value problem 

u rr + -u r + \u ee = (0 < r < 1, < 9 < 2tt) 
with a specified value at the boundary circle 

u{i,e) = f{9). 



Practice: What can we try? We have traveling waves, product solutions, 
power series. Anything else? 

Lets try separation u(r, 9) = R(r)F(9). You get 

R"F + ^-R'F + RF" = 0. 

We need to divide by F at least: 

1 1 F" 
R" + -R! + -^R— = 0. 

2 

That is not quite separated. Multiply by ^: 

r 2 R" + rR! F^_ Q 
R + F 

Now what do we have? The R part doesn't look at all familiar. How about 
the Fl Well, there is no doubt about that. We need F" to be a constant 
multiple of F, and we know that equation very well by now. The possibilities 
for F(9) are cos(a0), sin(a6>), e a6 , e~ ae , and a9 + b, for various numbers a, 
b, or certain linear combinations of those. Since we want u continuous in 9 
if possible, it has to connect smoothly when 9 loops around to 2ir which is 
the same angle as 0. So we better not use a9 + b unless a = 0. A constant is 
ok, say b = 1 for definiteness. Also any of the sin(n#) or cos(n#) for n > 1 
would be alright. So far our list of acceptable functions F is 

1, cos((9), sin((9), cos(26>), sin(26>), cos(30), sin(3(9), ■ ■ ■ 

and these will be enough for our purposes. We'll abbreviate the list as 
cos(n#), sin(n#), (n > 0). Then ^p- = — n 2 . 
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Practice: 1. Convince yourself of that — n 2 . 

2. Would F{9) = sin(i#) be continuous at if we are thinking of running 6 
from to 27r? 




- n 2 = 0, 



r 2 R" + rR' - n 2 R = (n > 0) 



Now what? This is a second order linear differential equation, but look at 
the coefficients. They are not constants. So we can't use the characteristic 
equation method; that is based on trying exponentials and it won't work 
with variable coefficients. 

Practice: Check that in Lecture 10 if you don't remember. 

We can still try something else. We haven't tried a power series for a while. 
As it happens, the solution for this is the easiest power series you've ever 
tried. It is left for the Problems. In the next lecture we'll continue with u. 

Problems 

1. Having no other good ideas about the unfamiliar R equation r 2 R" + rR' — n 2 R = 0, 
we could try a power series 

R(r) = a + air + a 2 r 2 H 

But wait. It works out so nicely I would like to save you some work. Try just 

R(r) = r q 

If that doesn't work, you can always try something harder. You will like this one. Trust 
me. 

2. Combine your solution for R in Problem 1 with its corresponding F(8) that we found, 
and check it in the Laplace equation to make sure that you really do have a solution. Try 
it for a couple of different values of n. 



Today: Putting together our Laplace solutions, we are lead to challenging 
questions. 

We are solving the boundary value problem 



39 Laplace leads to Fourier 



l l 

U rr + -U r + ^UQff 



= 



(0 < r < 1, < 9 < 2vr) 
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with a specified value at the boundary circle 

u (i,e) = f(9) 

and so far we have a list of solutions to the Laplace equation: 

r n cos(n6), r n sm(n6), (n > 0) and also 1. 
Now we have to think about the boundary conditions. 

Practice: What does the boundary look like for the case u(r, 6) = r sin(#)? 
What is rsin(#) anyway? 

Example: Suppose the boundary condition is that 
u(l,6>) = f(6) = 6cos(36>) 
Of our list of solutions to Laplace, the one which stands out is 

r 3 cos(36>) 

It is almost what we need. Would it be alright to multiply it by 6? (yes, see 
Problem 1.) So one solution to the BVP is 

u(r, 6) = 6r 3 cos(3<9) 

Some people call the graph in Figure 38.1 a saddle, and this r 3 cos(36>) a 
monkey saddle. Can you see why? 

Example: If the boundary condition is now 

u(l,0) = f{0) = 6cos(3<9) + 3sin(66>) 

what will the solution be? We can try a linear combination 

u(r, 6) — cir 3 008(36*) + C2r 6 sin(6^) 

At least that will be a solution to the Laplace equation (see Problem 2). 
Then we have to choose the coefficients Ck if possible to get the boundary 
values. It will work to take ci = 6 and C2 = 3. Answer: 

u(r, 0) = 6r 3 cos(3<9) + 3r 6 sin(60) 

Look, you just throw the r n factors in. After all this buildup, does that seem 
strange? 
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Fourier 

M. Fourier observed the examples we've seen. He felt that something is 
missing. Of course you probably feel too that our boundary conditions are 
somewhat artificial, made to fit the product solutions we found. Fourier 
wondered if you could work with some such boundary conditions as 

u (i,e) = f(e) = e 2 

or something like that, which has nothing to do with the cosines and sines. 
It is an interesting question to ask. Could it be possible that you could 
somehow expand / in terms of cosines and sines? 

Could a function f{6) which is not expressed in terms of the cos(n0) and 
sin(n0) actually have those hidden within it? 

The question is more strange too, because it is nearly the opposite of what 
we are used to, like 



cos(0) = 1 - 




That goes the wrong way-we wanted 2 in terms of trig functions. Eventu- 
ally we will achieve this on page 168. But if the answer to the question is 
yes, then you can solve the membrane problem by just putting in the powers 
of r as we did above: 

If f{Q) = c + 01 cos(0) + 61 sin(0) + a 2 cos(20) + b 2 sin(20) + ■ • ■ 

then you can solve the BVP as follows, assuming convergence: 

u(r,8) = c + r(ai cos(0) + 6isin(0)) + r 2 (a 2 cos(20) + 6 2 sin(20)) H 

Fourier wanted to know how to extract the frequencies hidden in all func- 
tions. 




Figure 39.1 What frequencies could be hidden in this function? It is plotted 
on the interval [0, 107r]. 
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Problems 

1. We need to check carefully that when u is a solution to the Laplace equation, then 
so is 6u, or cu for any constant c. How can that be checked? Start with the first order 
derivative 

du 
dr 

Is it true that = c^? Then how about the second derivatives, say (cu)ee7 Is that 
the same as c(ugg)? That is the essential idea behind the fact that 

V 2 (ctt) = cV 2 u 

Why does that prove what we need? 

2. We also need to know that a sum of harmonic functions is harmonic. Reason that out 
similarly to Problem 1. 

3. We also need to understand that a linear combination of harmonic functions is harmonic. 
Part of the calculation goes like this. 

V 2 (ci1ii + C2U2) = V 2 (ciUi) + V 2 (C2M2) = CiV 2 (lil) + C2V 2 («2) 

Which part of that is by Problem 1? 2? Why does that prove what we need? 

40 Fourier's Dilemma 

Today: We look into pictures to help Fourier find his hidden frequencies. 

We now know that it is easy to solve the boundary value problem 

u„ + - u r + \ u e e = (0 < r < 1, < 6 < 2ir) 

with a specified value at the boundary circle u(l,9) = f(0) provided that we 
can express f(9) in terms of a series made from our list 

cos(n#), sin(n0), (n > 0) and also 1. 

Fourier wanted to know how to extract the frequencies hidden in /. 




Figure 40.1 What frequencies could be hidden in this function on [0, 107r]? 
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Figure 40.2 How about now? The function has been modified a little to 
make some features stand out. Compare with the next figure. 




Figure 40.3 Two plain sinusoids, cos(#) and cos(49). The graphs are drawn 
together, but the functions are not added. One is nearly right and the other 
is perfect. Which is which? 

You see from the figures that you can sometimes visually estimate the hidden 
information if it is not too deeply hidden. Fourier wanted to dig further. 
We'll dig all the way to the bottom in the next lecture. For now we take 
special identities that you know about. 

Example: 1. You probably remember cos 2 (0) = \ - \ cos(20). 

2. From the complex exponential identity e 3lS = (e lS ) 3 we can derive that 

cos 3 (<9) = ? cos(0) + * cos(36») 

So there are at least two functions not originally given as a linear combination 
of the cos(n^) and sin(nS), for which we have a Fourier series. There are 
just two terms in each. 

We can think of the numbers 1 and 3 as the hidden (circular) frequencies 
and the coefficients | and ^ are hidden amplitudes in the function cos 3 (#). 
We got those examples sort of by luck, by trig identities. In the next lecture 
we will find out a systematic method for getting this information. 

Problems 

1. Try to answer the captions in the figures. 

2. Sketch a graph of the cos 3 ((?) and sketch the two terms in its Fourier series separately 
to see how they fit together. A linear combination of sinusoids is not usually a sinusoid. 
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41 Fourier answered by Orthogonality 



Today: Finally: How to find the hidden frequencies. 

We have seen only two cases where a function which was not originally 
expressed as a linear combination of sinusoids in fact is one. One example 
is 

cos 3 (fl) = -cos(0) + Jcos(3fl) 

But we also know there are good reasons to find these hidden frequencies 
if they exist, to solve the Laplace equation for example. If there were only 
two examples, we wouldn't bother. 

After Fourier's time a concept was developed to extract the hidden informa- 
tion from other functions. It is called "orthogonality" of functions. It works 
like this. People knew that 




cos(#) sin(0) dO = 



because an antiderivative is |sin 2 (#). 

Practice: Check that. 

And there are a lot of other similar integrals which are 0. 

Example: You don't really need an antiderivative, because cos(#) sin(d) is 
an odd function: the graph has the kind of symmetry where each positive 
part to the right of is balanced by a negative part to the left of 0, and vice 
versa. So it must integrate to 0. By the same argument, 

/TT 
cos(n6>) sin(m0) d9 = 
-TT 

for all n and m. 
Functions / and g are called "orthogonal" if 

^ f(9)g(9)de = 

J — TT 

You can compare this to the dot product of vectors. a\i + a 2 j + a%k and 
b\i + \>2 J + bsk are orthogonal if 

aih + a 2 b 2 + a^bs = 
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The integral is a sum, or at least a limit of sums, and the values of the 
functions sort of play the role of coordinates, so there is the comparison. 

This was a very new idea, to use geometric language to discuss functions. 

Example: We saw above that all the cos(n#) are orthogonal to all the 
sin(m#). It is also true that all the cosines are orthogonal to each other. 
That is harder to prove, but can be done using trig identities. Also all the 
sines are orthogonal to each other, and all those trig functions are orthogonal 
to the constant function 1. 

So there is a long list of functions all orthogonal to each other, 

1 , cos(fl) , sin(0) , cos(2#) , sin(2#) , cos(3#) , sin(30) , cos(4#) , sin(4#) 

sort of like i, j, k in 3-dimensional space. We have a lot more than 3 
orthogonal functions though, since in effect we are now working in an infinite 
number of dimensions. This is a new concept, part of our new language. 

But what else does that get you? It gets you everything. Why: Suppose it 
is possible to write a function / as 

f(Q) = c + ai cos(0) + bi sm(6) + a 2 cos(26») + b 2 sin(20) H 

but that we do not yet know what the hidden coefficients Cn, a n , b n are. 
Here is a way to find them. Think about the vector case, where a vector v 
can be written 

v = a\T+ a 2 j + a^k 

but somehow we do not know what the coefficients are yet. Dot with j, for 
example. You find the number a 2 by doing a dot product: 

v ■ j = ait- j + a2j- J + a^k - J=0 + a2 + = a2 

by the orthogonality and the fact that J* is a unit vector. So we can try the 
same thing with functions, for example 

f{0) cos(20) d6 = 

-7T 

c cos(2<9) dO + / ai cos(6>) cos(26») dO + / bi sin(6>) cos(26>) dQ 

-7T J — TT J —TT 

/7T PIT PIT 

a 2 cos(26>) cos(26>) d6+ j b 2 sm(26) cos(26») d6+ / a 3 cos(36) cos(26l) d0+- ■ ■ 
-TT J — TT J — TT 

/TT 
cos(26») cos(26») dO + + H 
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We still need to know that remaining integral cos (26) cos (20) dO. Its 
value is n, as can be worked out using trig identities. 

These ideas lead to the Fourier coefficient formulas 

a n = - [ f(6) cos(n#) d9, 



and 



b n = - r f(0) sin(nfl) dO, 

TT J- n 



Of course we only derived the 02 case but the ideas work for all. We have 
shown that if / has a Fourier series then these must be the coefficients. It 
is a uniqueness argument. 

It is for another course to study which functions actually have Fourier series. 

Problems 

1. Try a short Fourier series such as 

f(0) = cos(0) + O.3sin(50) 

Sketch the graph of the two terms and of / to see how they fit together and how the sum 
no longer has a graph like that of a trig function. 

2. Same as Problem 1 for 

f(0) = cos(96») + cos(lO0) 

This one is hard to do by hand, but is interesting because you get to see an example of 
amplitude modulation, AM. 

3. If you are a musician you might know that cos(2i) sounds an octave higher than cos(t), 
and cos(3t) is a fifth above that. Try to solve the wave equation for a vibrating string 

ytt = yxx 

y(0,t) = 

y(Tr,t) = 

y(x,0) = f{x) 



supposing the initial string shape f(x) has a Fourier series consisting only of the sin(n:r) 
terms. You may use the product solutions cos(nt) sin(na;) that we found in Lecture 30 for 
the wave equation. 

4. Work out the Fourier coefficient integrals for the function f(0) = 2 . You'll need to 
integrate by parts twice. The result ought to be the series 

9 2 = y - 4cos(0) + ^ cos (20) - A C os(3<9) + ^ cos(40) - ^ cos(5^) + ■ ■ ■ 
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It is for another course to discuss the convergence of Fourier series, but this one is alright 
for — 7T < 6 < 7r. Of course this series cannot be correct for all values of 8 because the 
series is periodic, and 8 2 is not. This series answers Fourier's question on page 162. 

5. In problem 4 evaluate the series at 9 = 7v to discover the value of Yl ^2 that we have 
needed ever since page 135. 

Can you find other interesting series? 

42 Uniqueness 

Today: How do we ever know that we have found the only solution to a 
problem? 

Some problems have only one solution. 

Example: If 3x - 2 = 5, then 3a; = 7, so x = |. That is a proof of 
uniqueness. 

Most people in reading that example think they have found the solution. 
But that is not quite true. We have to do this: 

Example: If x = |, then 3a; = 7, so 3x - 2 = 5. 

The problem is so easy that we automatically know all the steps are re- 
versible, and we feel no need to check it. 

Example: If there is an identity of the form 

\x\ = ax + (x- l) 4 + (x - l) 6 + (x - l) 8 
then take x = 1. You get 

a = 1 

That is a proof of uniqueness. 

Hmm. As far as that goes, it is ok, but something is badly wrong isn't it? 
There is no such identity of course, as you see by trying x = 0. 

Other examples: 

• We have a uniqueness theorem on page 21 about some ordinary differ- 
ential equations. That is for another course. 
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• On page 39 we proved by an energy argument that our favorite second 
order equation y" = —y has only those solutions that we know from 
calculus. 

• There is no choice about Fourier coefficients as shown in Lecture 41. 
If a function has a Fourier series, there is only one. But we avoid the 
question of when a function has a series. Too advanced. 

a Wave Equation example 

Here is a case where you can gain a good tool for the toolbox. This is an outline of an 
elegant argument which is usually considered to be outside the scope of this course. But 
some of you asked about it, and it is a good challenge. This is the wave equation for a 
string of length L 

Utt = tin 

u(0,t) = f(t) 
u{L,t) = g(t) 
u(x,0) = h(x) 
u t (x,0) = k(x) 

The initial shape h and velocity k are both given, and you can even shake both ends of the 
string by / and g. I claim that uniqueness holds: if U\{x, t) and u 2 (x,t) are solutions, then 
they are identically equal. How can this be proved? Let w be the difference, w — 112 — Mi . 
You'll need to convince yourself that subtracting wipes out all the boundary and initial 
conditions so that w solves the problem 

w u = w xx 

with 

w(0, t) = w(L, t) = w{x, 0) = w t (x, 0) = 
Once we prove that w — then we will have u\ = U2 like we claim. 

New idea: Let E — f () L + dx which we call the energy of the string. In a 

nutshell here is what happens: 

i) The energy is conserved i.e. constant in time. 

ii) E starts out being at t — 0. 

iii) The only way E can be identically is if w is identically 0. 

Of those items only i) needs a calculation, and the others just need some careful thinking. 
Here is the calculation without any explanatory comments: 

L r-L 

(w t wu + w x w xt ) dx = (w t w xx + w x w xt ) dx 
Jo 

L 

{w t w x ) x dx = [w t w x ] x=0 = 
Problem Fill in as many details of this argument as you can. 



dE_f 
dt ~ /„ 
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43 selected answers and hints 



page 4 

3. A function increases where the derivative is positive, so y is increasing on (—1,0) and 
(l,oo). 

4. If y(0) = 1/2 then y decreases as long as it is positive, so probably y — > 0. 

5. y' = (.03 + .005% 
page 7 

1. — kx~ 2 — mi", the minus because the force acts toward the large mass at the origin, 
the k as a generic constant, and m the small mass. Or, if you know the physics better, 
k — mMG where M is the large mass and G the gravity constant. G is about 10~ 10 [mks]. 

page 11 

2. On a left face of one of these, y x > 0. The equation tells you then y t > 0. So these 
waves must move left. 

3. The idea here is to find the derivatives. This is not strictly necessary because we 
already proved that waves traveling right at speed 3 are solutions no matter what shape 
they are. But it helps to check specific cases when you are learning something new. 

5. You can tell the velocity of the (model of the) dune, because that is how we interpret 
our solutions. But looking at the derivation of the model, we never connected the wind 
speed specifically with the dune speed. So you can't tell how fast the wind is. It is not 
part of the model. 

Naturally the real wind is going the same way as the real dunes, but not at the same 
speed. 

page 14 

1. Yes, because z'(t) = y'(t — 3) by the chain rule, and because y(a — y) evaluated at t is 
equal to z(a — z) evaluated at t— 3. The point here is that no explicit t occurs in the logistic 
equation. This time translation does not occur with an equation such as x' — xt — x 2 . 

3. Here is one way to outline the arithmetic. Write x = e~ atl as suggested. Starting with 
y = a(l + ae-^y 1 you have in years 1800, 1820, and 1840: 

5.3(1 + ci) = a 
9.6(1 + ax) = a 
17(1 + cix 2 ) = a 

The first and second expressions are both a, so they are equal, and this will give you x 
in terms of ci . Similarly the first and third give you x 2 in terms of ci . Consequently c\ 
solves a quadratic equation, which you can solve. Then it is easy to get a, x, and t\ in 
that order. 

page 17 

2. x(t) = -l/y/t + T/9 



170 



3. x(t) = l/y/t + 1/9 

10. x(t) = x(0)eti a{H)dt ^ 
page 20 

2. Trying u(x, t) — f(x)g(t) in u t + u x + u — gives fg' + f'g + fg = 0. It will separate 
if you divide by fg: 

I rt 

- + ^ + 1 = 

g f 

Since ^-M depends only on t, but it is equal to — y^y — 1 which does not depend on t, it 
must be a constant. Call it k. Then g' — kg and /' = — (k + l)f. Now apply what you 
know from problem 1. The waves don't all travel at the same speed. 

3. It is important to see cases where a method does not work, in addition to cases where 
it does. Proceeding as in problem 2 you get 

9 f 

which is not separated. There doesn't seem to be any way to separate this one. So try 
something else. This is normal, trying something else. 

4. We already know what happens to u t + cu x when u is a wave traveling right at speed 
c. So that isn't possible unless u = 0. But if you try u(x, t) = f(x — at) you find that you 
must have (c — a)f' + f 2 = 0, a separable equation. One of the traveling wave solutions 
is —l/(x — (c + l)t). Actually there are waves traveling right at every speed except c. 

page 24 

1. Once you get octave downloaded or matlab bought, and installed, start it. Once you're 
in, type 'help' with no quotes, press the enter key, and follow the information given. There 
is also 'help demo', 'help plot', and many other things which the first 'help' will tell you 
about. It is important to keep trying. At first, just type simple things like '2+3' or 'cos(pi)' 
until you get comfortable with it. After that, learn how to make an "m-file": that is a 
file called something like 'myfile.m' which contains just the text of some commands. For 
example, df ield is a large m-file for matlab, which doesn't seem to work in octave (but 
there is the df ield applet, too). Here is a much shorter example. 

2+3 
t=0: .1:3.5 
y=cos(t) ; 
plot(t,y) 

Then use the file by typing 'myfile' at the prompt. All the commands will be executed. If 
an error occurs just edit your m-file. This is easier usually than fixing any errors at the 
prompt. Once you can do this successfully, you are over the hurdle and you can use the 
software to help you learn. It is a nice tool to have, once you get past the first stage. The 
freeware xpp, or xppaut, is also excellent. 

page 28 
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1. This is tough because you have to integrate (e 1 ' 3t x)' — e 13t cost. The answer is 
x = ce 1 ' 3 ' + (-1.3cost + sint)/(l + (1.3) 2 ). 

5. x' = .03a; + 3000, x(Q) = 0. Multiplying by integrating factor e~ 03t you find 
(e- 03t x)' = 3000e- 03 '. So e - mt x = -lOOOOOe^ 03 * + c, x = -fOOOOO + ce 03t . Then 
x(0) = -100 000 + c = so c = f 00 000 and x(t) = 100 000(-1 + e 03t ). Finally, 
x(20) = 100 000(-l + e' 6 ) = 82 212. Not bad. 

6. Here is one philosophy. We have x' = —.028a; + .028 where x(t) is the balance on an 
account earning 2.8% interest and deposits of $.028 per year, and y = .028(y — y 2 ) where 
y(t) is a population with a natural growth rate of 2.8% and a limiting size of 1 (million 
or whatever). Change the units of x to dollars per person with the idea that the analysis 
certainly could be applied to one individual at a time. Change the units of y to people 
per dollar to acknowledge the fact that the economy determines how many people can 
be supported. Then the substitution x = 1/y is suggested by the units. This argument 
doesn't prove that a solution to one equation goes over into a solution to the other! Only 
the substitution can show whether that is true: if y = .028(y — y 2 ) and x = 1/y then 
x' = -y~ 2 y' = -y- 2 (.028(y - y 2 ))= -.028(y- 1 - 1)= -.028(z - 1). 

page 34 

3. sec 2 y dy — dt, so y = tan t. 

6. Yes. To get rid of the cosine you would have to start with a differential equation for 
tan -1 which contains only arithmetic operations H — */. We know that d(tan -1 t)/cK = 
1/(1 + t 2 ) so all you have to do is solve that numerically. 

8. a ■ 10 = b ■ 10 4 . So a = 10006. Then Newton says u' = -WOObu and Stefan says 
u = —bu 4 . The ratio is 1000u~ 3 so if u < 10 Newton cools faster. If u > 10 Stefan cools 
faster. 

page 39 

1. s" + s 2 = (xi + x 2 )" + (xi + x 2 ) 2 = x'i + x'i + x\ + 2x!X 2 + x\ = + + 2a;ia;2. Uh-oh. 

3. No, s solves s" + s = 2cos8t instead. 

8. x'(x" + x 3 ) = integrates to + ^ = c. 

10. Somebody tried to use the characteristic equation for a nonlinear equation, which is 
very rong. Then even if r 2 + 1 = had been correct, the roots should have been ±i. 
Finally, you should never just add "+c" to something at random. 

page 42 

1. T(t) = cie~ 2 \ X(x) = c 2 cos(^a;) + c 3 sin(V2a;), u(x,t) = e~ 2t (c 4 cos(^a-) + 
C5 sin(v / 2a*). 

3. You have T" = which gives the steady state solutions, 
page 47 

2. Trying y = e rt we need r 2 + 3r + 4 = 0, r = (-3 ± y/9 - 16)/2 = -3/2 ± So 
y(t) = e- 3t/2 (a cos(\/7t) + c 2 sin(\/7t)) 

3. If a + i = then a 2 + 1 = 0. So a = ±i. Checking, those both work, since \ = — i. 
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If a + \ = i then a 2 — ia + 1 = 0. Quadratic formula gives a — 2 1 4 . 

6. There is a factor of e J so this grows, and a factor of cos(t) + i sin(i) which goes around 
the unit circle counterclockwise, so this is a spiral. 

7. r 2 + br + c — (r + 2 + i)(r + 2 — i), so multiply it out. 
page 48 

4. If you aren't sure, try some examples, but the identity 

p(x) = do + aix H + ctn-ia;™ -1 + x" = (x - n)(x - r 2 ) ■ ■ ■ (x - r„) 

gives n + r 2 H + r„ = -o„_i. 

page 53 

1. Csin(/t + 0) = C sin(ft) cos(0) + C cos(ft) sin(0). So to match this with the given 
A cos(/i) + B sin(/t) you have to solve C cos(<f>) = _B and C sin(0) = A for C and 0. That 
is like saying that (B, A) is a point in the plane and you want to find polar coordinates 
(r, 6) = (C, 0) for this point. That is always possible. 



2. This is a little harder, because the only clue you have is that in the given identity, the 
fY-o,-.i,or,^;c, s |ij> anc j ~L are half the sum and difference of the frequencies ^ and We 



have 

sin(a + 6) = sin(a) cos(fo) + cos(a) sin(b) and sin(a — b) = sin(a) cos(6) — cos(a) sin(b) 
so 

A sin(a + b) + B sin(a — b) = (A + B) sin(a) cos(b) + (A - B) cos(a) sin(6) 
For example sin(8t) + 2.3 sin(7t) = 3.3 sin(7.5t) cos(.5t) - 1.3cos(7.5t) sin(.5t). 

3. From problem 2 we see that when you graph something like sin(7.5t) cos(.5t) you see 
a slow and a fast frequency in spite of the fact that you have started with the two near 
frequencies ^ and ^ . In the case of the Figure the near frequencies are ^ and ^ , so 
the corresponding slow and fast frequencies are 

1 /\/2 

2 V 2tt 2tt 

which are about .12 and .35. The situation is more complicated than the identities in 
problem 2 though, because you really need the awful identity A sin(a+b)+B sin(a— b+(j)) = 
(A+B cos(0)) sin(a) cos(6)+ (A— B cos(0)) cos(a) sin(&)+ (B sin(0)) sin(6)(cos(a) — sin(a)). 
The only important thing is that the slow and fast frequencies still show up. 

page 57 

1. X2 — xi — x[ from the first equation, then the second gives x[ — x'{ = x\ + 2(xi — x[), 
or x'{ — 3x[ + 3xi = 0. 

2. It is the same as for x\. Either calculate it, or recognise that the formula X2 = x\ — x\ 
expresses x 2 as a linear combination of two solutions to the equation for x\ . 

4. Assuming y" — 3y' + 3y = we have zi = y' — y = z[ — z\, z 2 = (y' — y)' = y" — y' = 
3y' — 3y = 3z 2 . So this system is 

z[ = Zl + Z 2 
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Z*2 — 3,?1 

This is different from the system of problem 1 even though they have the same second- 
order equation associated with them. The point is that there is nothing unique about a 
system associated to a higher order equation. 

page 61 

1. If you use y = x' then y = x" = x 3 — x so the system is 

„3 



y, 



y 



But if you take say y — x' + x, then y' — x" + x — x 3 — x + (y — x) — x 3 — 2x + y 

6. A conserved quantity is 2(x') 2 + x 4 . The level sets are noncircles whose shape changes 
depending on the value. 



page 66 

2. A(u + v) = Au + Av — b + b. So a sum of solutions to the equation Ax 
usually a solution. Usually 26 / b. 



b is not 



But if Av = you get A(u + v) = Au + Av = b + = b. So adding a solution of Ax — 
to a solution of Ax = b gives a solution to Ax — b. 

And, as a special case, adding two solutions of Ax — gives a solution to Ax = 0. 

3. All the drinks must go to New York or Chicago, so 9000 = 200y + 180a;. All the fuel 
must be used on the way, so 260000 = 6500y + 5400a;. These are our system. It seems 
fairly clear that this is too simplified to be realistic. If we really wanted to model costs of 
an airline we would certainly not consider drinks as important as fuel. The next largest 
expenses after fuel would have to be identified. They might be the cost of the planes, or 
maintenance, or personel salaries and fringes. 

4. This kind of product has not been defined, 
page 70 

1. row 1 has been replaced by (row 1-1.6 row 2). 



3. The reduced form is 



5.5 
.5 



4. The third equation says = 2 which is not true. So some incorrect assumption must 
have been made when the equations were first used. Be sure you get this: When we 
started doing row operations, we assumed there was a solution. So there are no solutions. 

8. That is one solution, not 'no solutions'. 

page 74 



3. 
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1 



with no restriction on a. For the 4 by 4, if a / then 



If a — then the 4 by 4 is not invertible. 
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For example, it can't be invertible because 



10 0] I'll TO' 

1 

1 





5. Ix — x for all x, but is not an eigenvector. So, all vectors except are eigenvectors 
of I, with eigenvalue 1. 

6. The octave command works correctly on the case which has a unique solution, and 
seems to give some answer for the other two cases. 

page 79 





X 




"1 2 


3 




"1 


2 


3" 


3. We have to solve A 


V 


= 0. Row operations give A — > 


-3 


-6 







1 


2 




z 




-6 


-12 














1 -1 
1 2 




so y — — 2z and x — z. So A 



0. 



8. For this problem and for problem 7, note that A and 7 'A have the same eigenvectors, 
but the eigenvalues of 7 A are 7 times as large. 

11. The calculation here is correct but misguided! As soon as you know that det(A — XI) — 
(3 — A) (7 — A), stop, because this is already factored. There is no need to multiply it out 
and use the quadratic formula. 



page 82 

1. This matrix gives y — Ax = 



8x4 

Xl 
X2 
X3 



and it is easy to invert that: x — 



J/2 
J/3 



But if 



you want to calculate the determinant as a check, the determinant is —8 / 0. 
page 86 

2. x(t) = 0, and y(t) = 5e 3 * 

5. x = ae 2t and y = fee -3 '. If some of a or b are you get solutions running along a half axis 
or stationary at the origin. Otherwise y~ 2 ^ = {const)e^ 2 ^ 3 ^^' i ^ = (const)e 2t = (const)x. 

9. x(0) =u + 2v, so x(t) = ue 5t + 2ve At . 

10. Don't add "+C" without a reason. This is correct otherwise, 
page 89 

3. The two plots are identical. 

5. The real and imaginery parts were taken incorrectly because somebody forgot that the 

vector \ is also complex. 
i\ 

page 95 
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3. All rational numbers have either terminating or repeating decimal expansions, while 
the irrational ones have the other kind. All integers are either positive, negative, or 0. All 
math teachers are mean. 

4. The system of equations is just x' = y' = 0, so the solutions are constant, and in the 
picture this means they don't go anywhere. 

5. Don't forget that the differential equation tells you a vector field, so it shows you which 
way the spiral goes. 

page 99 

3. I tried the system x' = x — xy, y' = —y + xy — .2. I found that solutions can make 
increasing oscillations until the prey becomes extinct. This is not what one would hope 
for! My advice is not to fool with Mother Nature. 

6. (0,0) only, 
page 102 

5. The correct conclusion is that chaos is possible in a system of two or more masses, 
page 106 

1. For c = there is the solution y = 0. For c < the general solution involves linear 
combinations of e ±v ^ x and these all get wiped out by the boundary conditions, except 
for y = 0. 

4. If c < then the only answer is y — 0. If c > write it as c = d 2 where you may as 
well assume d > 0, and the general solution to y" = —d 2 y is y(x) = acos(dx) + bsin(dx). 
The boundary condition y(0) = forces a = 0. Then the condition y'(3) = = db cos(3d) 
forces 3d = f , 'f, ^f, . . . . Thus y{x) = 6cos (!f), for k = 1, 3, 5, . . . . No restriction on b. 

page 111 

2. In Mt = au xx — u, the — u term certainly seems to be a negative influence on the time rate 
of change. So it seems fairly reasonable that Mars bars will cool faster than Earth bars. 
That argument doesn't work if u is ever negative. Suppose that u is 'bounded below', i.e. 
u > — m for some number m. Then replace tibyii+mto get positive temperatures only. 
That fixes the argument. So Mars bars cool faster than earth bars unless they have some 
with infinitely negative cold temperatures. 

page 115 

1. This problem concerns heat conduction along a bar which is initially 400 degrees, and 
which has insulation all around, including the ends. The steady state temperature is what 
you think it is. 

5. 3cos(2;)e~* 

10. the heat equation for u = e~ at y says that u t = —ae~ at y(x) = u xx = e~ at y"(x). So 
we must have 

y" = -ay with i/(0) = y(ir) = 

page 120 

3. u(x, t) = e~* sin(a;) + .5e~ 32 * sin(3t) + .25e -52 * sin(5i) 
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page 125 

3. sin(x + 3i) + sin(:r — 3i) = sin(a;) cos(3i) + cos(:c) sin(3i) +sin(a;) cos(3i) — cos(:r) sin(3i) = 
2 sin(:r) cos(3t) has the initial condition 2sin(:r) of problem 2. For this problem the only 
difference is that we want 2sin(2a;) instead, so the solution is 2 sin(2x) cos((2)(3i)). You 
have to multiply the 3t by 2 because otherwise you would not be creating traveling waves 
which are functions of x ± 3t, as required by the wave equation utt = 3 u x 



' XX ■ 



8. We have said only that sums of left and right traveling waves are solutions. In this 
case the function is seen to be a solution to utt = u xx by differentiating it, or by using the 
idea of problem 3 to write it as a sum of left and right traveling waves: cos(t) sin(x) = 
(l/2)(sin(x + t) + sin(x - t)). 

page 129 

3. cos(fcn 2 i) sin(na:) where k = \J ^y, n = 1,2, 3... 

4. the answer is not unique because of insufficient BC 
page 135 

2. x = \ works, so s is no more than 2. 
page 140 

3. The radius of convergence is no affected by the —10, so it is 1, just like in Problem 
1. But when x = 8.3, \x — lOj > 1 so you are outside the radius of convergence. You are 
trying to add powers of 1.7, which diverges. 

page 148 

1. Jo (4) = s 6 (4) + E = l-f+E = -l+E 

2. The tail E = 2 2 4 2 6 2 g 2 — 2 2 4 2 6 2 s 2 io 2 + " ' = 2 2 4 2 6 2 s 2 [} — To 2 + io 2 i2 2 — " ' ) ^ y° u 
remember about alternating series, that is less than 2 a 4 2 t 6 2 8 2 — §• 

So Jo (4) < -| + | < as needed, 
page 157 

1. Integrating 717 + 287 % = dS you find In (r 71 V 287 ) = S + c 
page 164 

1. In the figure there seem to be at least two frequencies at work. The larger one has 10 
cycles in the interval [0, Wtt] and seems to start positively at 9 — 0. So a multiple of cos(#) 
is a good guess for that one. The smaller one is harder to count. It has less amplitude 
and about 4 cycles for each of the larger ones. It is probably cos or sin of 49 or 59, with 
a smaller amplitude than the other one. 

page 167 

3. We know about the product solutions cos(ni) sin(na;) for the wave equation with these 
string boundary conditions. And we know that you can form linear combinations of these 
solutions. So the idea: If we could write / as a series of sines, then we could use those 



177 



same coefficients in forming our linear combinations. If 

f(x) — bi sin(x) + &2 sin(2a;) + 63 sin(3a;) + • • • 

try 

y(x, t) — b\ cos(t) sin(a;) + &2 cos(2t) sin(2:r) + 63 cos(3t) sin(3x) + • 
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